Luận án Tiến sĩ Cơ học kỹ thuật: Phân tích đáp ứng của profile cánh máy bay theo cách tiếp cận đối ngẫu

ĐẠI HỌC QUỐC GIA HÀ NỘI

TRƯỜNG ĐẠI HỌC CÔNG NGHỆ

--o0o--

NGUYỄN MINH TRIẾT

PHÂN TÍCH ĐÁP ỨNG

CỦA PROFILE CÁNH MÁY BAY

THEO CÁCH TIẾP CẬN ĐỐI NGẪU

Chuyên ngành: Cơ kỹ thuật

Mã số: 62 52 01 01

LUẬN ÁN TIẾN SĨ CƠ HỌC KỸ THUẬT

NGƯỜI HƯỚNG DẪN KHOA HỌC

1. GS.TSKH Nguyễn Đông Anh

2. PGS.TS Phạm Mạnh Thắng

HÀ NỘI - 2017

II

Tôi  xin  bày  tỏ  lời  cảm  ơn  sâu  sắc  tới  các  thầy  hướng  dẫn  khoa  học  là  GS.TSKH. Nguyễn Đông Anh và PGS.TS. Phạm Mạnh Thắng, các thầy đã trực tiếp  hướng dẫn tận tình và giúp tôi hoàn thành luận án này.

Tôi cũng chân thành cảm ơn các nhà khoa học và các cán bộ của khoa Cơ học kỹ  thuật & Tự động hóa, trường Đại học Công nghệ, Đại học Quốc gia Hà Nội, và Viện  Cơ học, Viện Hàn lâm Khoa học & Công nghệ Việt Nam, đã tạo điều kiện thuận lợi,  giúp đỡ tôi trong quá trình học tập, nghiên cứu tại đây.

Hà Nội, ngày … tháng … năm 2017

Nguyễn Minh Triết

LỜI CẢM ƠN

III

Tôi  xin  cam đoan đây  là  công  trình nghiên  cứu  của riêng tôi.  Các  số liệu,  kết

quả nêu trong luận án là trung thực và chưa từng được ai công bố trong bất kỳ công

trình nào khác.

Hà Nội, ngày tháng năm 2017

Tác giả luận án

LỜI CAM ĐOAN

Nguyễn Minh Triết

IV

LỜI CẢM ƠN ............................................................................................................ II   LỜI CAM ĐOAN .....................................................................................................  III   MỤC LỤC ................................................................................................................  IV   DANH MỤC MỘT SỐ KÝ HIỆU VÀ CHỮ VIẾT TẮT ........................................... VI   DANH MỤC HÌNH VẼ, ĐỒ THỊ  .........................................................................  VIII   DANH MỤC BẢNG VÀ SƠ ĐỒ KHỐI .................................................................... IX  MỞ ĐẦU ..................................................................................................................... 1  1. Tính cấp thiết của đề tài ....................................................................................... 1  2. Mục tiêu nghiên cứu ............................................................................................ 2  3. Đối tượng nghiên cứu .......................................................................................... 3  4. Nội dung nghiên cứu ........................................................................................... 3  4.1. Phương pháp nghiên cứu .............................................................................. 3  4.2. Hướng giải quyết .......................................................................................... 3  4.3. Kết quả dự kiến ............................................................................................ 3  5. Cấu trúc của luận án ............................................................................................ 4  CHƯƠNG 1. TỔNG QUAN VỀ BÀI TOÁN PHÂN TÍCH ĐÁP ỨNG CỦA THIẾT  DIỆN CÁNH ............................................................................................................... 7  1.1. Khái niệm cơ bản về khí đàn hồi ....................................................................... 7  1.2. Các nghiên cứu đáp ứng của thiết diện cánh ..................................................... 8  1.3. Thiết diện cánh phi tuyến ................................................................................ 13  1.4. Một số nghiên cứu liên quan ở trong nước ...................................................... 16  1.5. Cách tiếp cận đối ngẫu .................................................................................... 18  1.6. Vấn đề nghiên cứu của luận án ....................................................................... 19  Kết luận chương 1...................................................................................................... 20  CHƯƠNG  2.  MÔ  HÌNH  CƠ  HỌC  CỦA  THIẾT  DIỆN  CÁNH  CHUYỂN  ĐỘNG  TRONG DÒNG KHÍ ................................................................................................. 21  2.1. Lực khí động dừng và tựa dừng ...................................................................... 21  2.1.1. Lực khí động dừng .................................................................................. 21  2.1.2. Lực khí động tựa dừng ............................................................................. 25  2.2. Phương trình chuyển động của thiết diện cánh ................................................ 28  2.3. Hiện tượng flutter ........................................................................................... 30  2.3.1. Hiện tượng mất ổn định 1 bậc tự do ......................................................... 30  2.3.2. Hiện tượng mất ổn định 2 bậc tự do ......................................................... 32  2.4. Tính toán vận tốc flutter trong hệ tuyến tính ................................................... 34  2.4.1. Hệ tự dao động tổng quát ......................................................................... 34  2.4.2. Thiết diện cánh 2 chiều có điều khiển PID ............................................... 36

MỤC LỤC

2.5. Tính toán thiết diện cánh bằng phương pháp CFD .......................................... 37  2.5.1. Mô phỏng khí động lực trên mô hình cánh máy bay ................................ 38  2.5.2. Tối ưu hình dạng khí động sử dụng phương pháp SQP ............................ 45  2.5.3. Mô phỏng CFD trên cánh máy bay với các góc tới lớn ............................ 53  Kết luận chương 2...................................................................................................... 61  CHƯƠNG 3. PHÁT TRIỂN KỸ THUẬT ĐỐI NGẪU CHO BÀI TOÁN DAO ĐỘNG  PHI TUYẾN .............................................................................................................. 62  3.1. Phương pháp tuyến tính hóa tương đương ....................................................... 62  3.1.1. Tiêu chuẩn tương đương kinh điển .......................................................... 63  3.1.2. Tiêu chuẩn sai số thế năng ....................................................................... 64  3.1.3 Tiêu chuẩn tương đương điều chỉnh .......................................................... 65  3.2 Tiêu chuẩn đối ngẫu có trọng số ...................................................................... 66  3.3. Những cải tiến của phương pháp đối ngẫu có trọng số .................................... 68  3.3.1. Cải tiến 1 ................................................................................................. 68  3.3.2. Cải tiến 2 ................................................................................................. 69  3.3.3. Cải tiến 3 ................................................................................................. 69  3.4 Áp dụng cho dao động tự do của hệ phi tuyến dạng Duffing bậc cao ............... 70  3.5. Áp dụng cho dao động ngẫu nhiên .................................................................. 73  CHƯƠNG  4.  ÁP  DỤNG  KỸ  THUẬT  TUYẾN  TÍNH  HÓA  ĐỐI  NGẪU  CHO  BÀI  TOÁN PHÂN TÍCH ĐÁP ỨNG PHI TUYẾN CỦA THIẾT DIỆN CÁNH ............... 76  4.1. Mô hình thiết diện cánh .................................................................................. 76  4.2. Phương trình xác định vận tốc tới hạn ............................................................. 79  4.3. Áp dụng kỹ thuật tuyến tính hóa đối ngẫu ....................................................... 81  4.4. Các ví dụ và tính toán bằng phương trình vi phân ........................................... 84  4.4.1. Số liệu đầu vào ........................................................................................ 84  4.4.2. Tìm vận tốc tới hạn bằng phương pháp số................................................ 87  4.5. Kết quả tính toán với ví dụ 1 ........................................................................... 89  4.6. Kết quả tính toán với ví dụ 2 ........................................................................... 90  4.7. Kết quả tính toán với ví dụ 3 ........................................................................... 92  4.8. Kết quả tính toán với ví dụ 4 ........................................................................... 94  4.9. Kết quả tính toán với ví dụ 5 ........................................................................... 97  Kết luận chương 4.................................................................................................... 100  KẾT LUẬN VÀ KIẾN NGHỊ .................................................................................. 102  TÀI LIỆU THAM KHẢO........................................................................................ 107  PHỤ LỤC ................................................................................................................ 116

V

VI

véc tơ hàm phi tuyến

A

Góc tới hoặc góc nâng (Angle of Attack)

AoA

véc tơ, hàm tuyến tính tương đương

B

hệ số tuyến tính hóa tương đương

b, k

hệ số chuẩn hóa

C

Các hệ số khí động (nâng, cản, mô men)

CL, CD, CM

hệ số độ cứng tuyến tính

c1,ktt,kα

hệ số độ cứng phi tuyến

c3, c5

Thiết kế với hỗ trợ của máy tính (Computer-Aided Design)

CAD

Động lực học chất lỏng tính toán (Computational fluid dynamics)

CFD

Mô hình cấu trúc tính toán (Computational Structural Model)

CSM

Lực cản (Drag)

D

tỉ số các hệ số tuyến tính hóa theo các tiêu chuẩn

dk(µ)

),

xxD t ( , 1

t D    2

12

hiệp phương sai

  e x x   ,

sai số phương trình

E{.}, <.>

kỳ vọng toán học

f(t), u(t)

kích động ngoài

hàm phân phối xác suất

F(x)

Tương tác dòng khí kết cấu (Fluid Structure Interaction)

FSI

 g x x    ,



hàm phi tuyến của dịch chuyển

h

hệ số cản tuyến tính

H x x    ( , )

hàm tổng năng lượng

Phương trình Karush-Kuhn-Tucker

KKT

ma trận hệ số khuếch tán

K(x,t)

Các hệ số bộ điều khiển PID

Kp, Ki, Kd

Lực nâng (Lift)

L

Dao động vòng giới hạn (Limit Cycle Oscillation)

LCO

DANH MỤC MỘT SỐ KÝ HIỆU VÀ CHỮ VIẾT TẮT

m

khối lượng

minS

giá trị cực tiểu của tiêu chuẩn tuyến tính hóa

trung bình xác suất

mx

p

trọng số

p(μ)

hàm trọng số

P{.}

xác suất của một sự kiện

Bộ vi tích phân tỉ lệ (Proportional Integral Derivative)

PID

PTTH

Phần tử hữu hạn (Finite Element)

q

Áp suất khí động lực

r

hệ số tương quan

hàm tương quan

R(t1,t2)

Hệ Reynolds Navier-Stokes (Reynolds Averaged Navier-Stokes)

RANS

Re

Số Reynolds

S

biểu thức tính diện tích

Lập trình toàn phương liên tiếp (Sequential Quadratic Programming)

SQP

hàm mật độ phổ

Sx(ω)

T

chu kỳ dao động

t

thời gian

TMD

Bộ hấp thụ dạng khối lượng (Tuned Mass Damper)

TTH

Tuyến tính hóa (Linearization)

UAV

Máy bay không người lái (Unmanned aerial vehicle)

v t x t u, ( ), ( )

vận tốc

x(t)

dịch chuyển

 x t

gia tốc

X, Y

biến ngẫu nhiên

α, β

các hệ số hằng, hoặc góc tới

δ(x)

hàm Delta Dirac

θ

góc giữa hai véc tơ

VII

hệ số trở về

λ

mức độ phụ thuộc tuyến tính

μ

Khối lượng riêng

mô men trung tâm

μn

mô men liên kết trung tâm

μnm

độ lệch chuẩn

σx

2

x

phương sai

độ trễ

τ

tần số kích động

ω

tần số dao động tự do

ω0

VIII

IX

Hình 1. Tam giác khí đàn hồi của Collar

Hình 2. Mô hình thiết diện cánh hai chiều, lực khí động quy về lực tập trung

Hình 3. Mô hình cánh theo tấm bị ngàm

Hình 4. Mô hình cánh phẳng, lực khí động tính bằng CFD

Hình 5. Mô hình kết hợp CFD-CSM

Hình 6. Mô hình thiết diện cánh 2 chiều có cánh nhỏ điều khiển

Hình 7. Phi tuyến bậc ba và phi tuyến khe hở tự do

Hình 8. Dao động vòng giới hạn điển hình

Hình 9. Một số thuật ngữ về cánh

Hình 10. Dòng dừng đi qua một thiết diện cánh 2 chiều

Hình 11: Đường dòng và các đường đẳng thế của một xoáy 2 chiều tại gốc tọa độ

Hình 12: Dòng không dừng đi qua thiết diện cánh

Hình 13. Mô hình cánh 2 bậc tự do

Hình 14: Mô hình lực chất lỏng tác động vào hệ 1 bậc tự do

Hình 15. Các hệ số khí động theo góc xung kích

Hình 16. Mô tả hiện tượng mất ổn định lên xuống –xoắn

Hình 17. Phác thảo thiết kế thiết diện cánh máy bay.

Hình 18. Mô hình cánh máy bay tạo ra trong ANSYS.

Hình 19. Chia lưới trong ANSYS

Hình 20. Biểu đồ phân bố áp suất

Hình 21. Biểu đồ phân bố vận tốc

Hình 22: Ứng suất tương đương trên cánh

Hình 23. Quan hệ giữa lực nâng với vận tốc tương đối

Hình 24. Quan hệ giữa lực cản với vận tốc tương đối

Hình 25. Hình dạng thiết diện cánh máy bay ban đầu và cánh tối ưu

Hình 26.  Hình ảnh phân bố áp suất tại Re=5.105 và a=5o trên thiết diện cánh Eppler

66 (hình trái) và thiết diện tối ưu theo SQP (hình phải).

DANH MỤC HÌNH VẼ, ĐỒ THỊ

Hình 27. Hình ảnh phân bố vận tốc Re=5.105 và a=5o trên thiết diện cánh Eppler 66

(hình trái) và thiết diện tối ưu theo SQP (hình phải)

Hình 28: Hình ảnh phân bố ứng suất trên cánh tại Re=5.105 và a=5o đối với thiết diện

cánh  Eppler 66 (hình trái) và thiết diện tối ưu theo SQP (hình phải)

Hình 29. Phác họa thiết diện đặc trưng của cánh máy bay

Hình 30. Hình học của hình dạng cánh

Hình 31. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh

tại góc tới là 0 độ

Hình 32. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh

tại góc tới là 10 độ

Hình 33. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh

tại góc tới là 18 độ

Hình 34. Đồ thị hệ số lực nâng theo góc tới

Hình 35. Đồ thị hệ số lực cản theo góc tới

Hình 36. Mô hình khí động với các bậc tự do tịnh tiến và xoắn.

Hình 37. Mô hình giản đồ của bố trí thí nghiệm

Hình 38. Mô hình giản đồ của số liệu đầu vào và đầu ra cho việc đo và điều khiển mô

hình thí nghiệm

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 (cắt từ bài báo Li vcs 2011)

Hình 41. Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 1

Hình 42. Kết quả so sánh đường cong tần số-vận tốc trong ví dụ 2

Hình 43. Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 3, Kp=-0.5

Hình 44. Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 3, Kp=0.5

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)

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)

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)

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)

X

XI

Bảng 1. Các đặc điểm của Hợp kim nhôm 7075-T6

Bảng 2. Các đặc tính khí động lực trên các thiết diện cánh tại Re=5.105 và a=5o

Bảng 3. Các giá trị lực nâng và lực cản tương ứng với các góc tới khác nhau

Bảng 4. So sánh các tần số tính toán với tần số chính xác

Bảng 5. So sánh các dịch chuyển bình phương trung bình

Bảng 6. Số liệu của ví dụ 1

Bảng 7. Số liệu của ví dụ 2

Bảng 8. Số liệu của ví dụ 3

Bảng 9. Số liệu của ví dụ 4

Bảng 10. Số liệu của ví dụ 5

Bảng 11. So sánh một số vận tốc flutter trong ví dụ 1

Bảng 12. So sánh một số vận tốc flutter trong ví dụ 2

Bảng 13. So sánh các vận tốc flutter trong ví dụ 1

Bảng 14. So sánh các vận tốc flutter trong ví dụ 2

Bảng 15. So sánh các vận tốc flutter trong ví dụ 3, Kp=-0.5

Bảng 16. So sánh các vận tốc flutter trong ví dụ 3, Kp=0.5

Bảng 17. So sánh các vận tốc flutter trong ví dụ 4, Ki=-2 (1/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)

Bảng 19. So sánh các vận tốc flutter 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)

DANH MỤC BẢNG VÀ SƠ ĐỒ KHỐI

MỞ ĐẦU

1. Tính cấp thiết của đề tài

Hiện  nay,  máy  bay  là  phương  tiện  không  thể  thiếu  được  trong  mỗi  quốc

gia, có vai trò đặc biệt quan trọng trong lĩnh vực dân sự và quân sự. Do vậy, mặc

dù đã được phát minh và đưa vào sử dụng từ 100 năm trước, các nghiên cứu về

máy bay vẫn được tiếp tục phát triển mạnh trên thế giới nhằm nâng cao độ ổn

định, an  toàn và  tốc độ cho  các chuyến  bay.  Khi  máy  bay chuyển  động  trong

dòng khí sẽ xuất hiện các hiệu ứng khí động học, trong đó dao động flutter của

thiết  diện  cánh  máy  bay  rất  được  quan  tâm.  Phân  tích  đáp  ứng  của  thiết  diện

cánh máy bay là một bài toán quan trọng phục vụ quá trình thiết kế, chế tạo, vận

hành và bảo dưỡng máy bay. Thiết diện cánh chuyển động trong dòng khí không

nén  được  thường được mô  hình  bằng  mô hình  cơ  học hai  chiều.  Phương  trình

chuyển  động  ứng  với  mô  hình  thường  là  hệ  tự  dao  động  và  có  tính  chất  phi

tuyến. Khi nghiên cứu hệ phi tuyến này, người ta quan sát thấy hiện tượng mà ở

đó có xuất hiện vòng giới hạn, các hiện tượng rẽ nhánh Hopf và hiện tượng mất

ổn định flutter. Vấn đề khoa học này đã thu hút nhiều nghiên cứu trong những

thập niên trở lại đây, nhất là nghiên cứu phục vụ nhu cầu chế tạo các loại máy

bay với nhiều tính năng, đảm bảo ổn định khi bay ở các độ cao, vận tốc và điều

kiện bay khác nhau. Các  phương trình chuyển động của  thiết diện cánh đều là

phương  pháp  đã  có  để  có  thể  thu  được  lời  giải  đạt  được  độ  chính  xác  mong

muốn.  Mới đây, một cách tiếp  cận mới cho bài toán phi tuyến về dao động và

điều khiển kết cấu được đề xuất. Cách tiếp cận mới được biết với tên gọi cách

tiếp  cận  đối ngẫu với quan  điểm  tạo ra một  sự hài  hòa trong  nghiên  cứu,  cho

phép phát hiện bản chất của vấn đề một cách đầy đủ hơn. Áp dụng cho lĩnh vực

tuyến  tính  hóa  (TTH)  tương  đương  đã  dẫn  đến  phương  pháp  cực  tiểu  bình

phương  đối  ngẫu.  Ban  đầu  phương  pháp  được  đề  xuất  trong  nghiên  cứu  dao

động ngẫu nhiên của các hệ phi tuyến với kích động ngoài ồn trắng. Kết quả thu

1

phương  trình  phi  tuyến  và  có  thể  phi  tuyến  mạnh,  do  vậy  phải  phát  triển  các

được chỉ ra rằng với các hệ phi tuyến mạnh, phương pháp cho kết quả khá tốt và

phù  hợp  với  kết  quả  thu  được  từ  nghiệm  chính  xác  trong  trường  hợp  hệ  phi

tuyến có được nghiệm chính xác, và kết quả thu được từ nghiệm mô phỏng số

trong trường hợp không tìm được nghiệm chính xác của hệ phi tuyến. Tính ưu

việt của nó sau đó còn được tìm thấy trong trong nghiên cứu hệ nhiều bậc tự do

chịu  kích  động  ngẫu  nhiên.  Ý  tưởng  của  phương  pháp  được  mở  rộng  sang

nghiên cứu điều khiển giảm dao động cho hệ TMD. Các kết quả thu được về đáp

ứng của hệ TMD cũng tốt hơn hẳn so với các kết quả đã có trước đây. Cách tiếp

cận đối ngẫu ở trên có tính linh hoạt và có thể áp dụng được cho nhiều lớp hệ

phi tuyến khác nhau. Đây cũng là chủ đề của luận án với mục đích nghiên cứu

phát triển và áp dụng cho bài toán ổn định thiết diện cánh máy bay. Nghiên cứu

nhằm tìm ra những nghiệm gần đúng của bài toán với sai số nghiệm nhỏ so với

các nghiệm mô phỏng số trong trường hợp hệ đang xét có tính phi tuyến, thậm

chí là phi tuyến mạnh. Hướng nghiên cứu này chưa từng được triển khai cho đến

thời điểm hiện nay. Việc triển khai nghiên cứu sẽ tạo ra khả năng thu được các

kết quả mới chính xác hơn các kết quả đã biết, mở đường cho một cách tiếp cận

mới trong nghiên cứu các kết cấu hàng không và vũ trụ.

Như  vậy qua phân tích và nghiên cứu các tài liệu khoa học công nghệ có

thể  khẳng  định  đề  tài  “PHÂN  TÍCH  ĐÁP  ỨNG  PHI  TUYẾN  CỦA  THIẾT

DIỆN CÁNH THEO CÁCH TIẾP CẬN ĐỐI NGẪU” của luận  án có tính cấp

thiết, ý nghĩa khoa học và thực tiễn.

2. Mục tiêu nghiên cứu

-  Phát  triển  phương  pháp  luận  cho  cách  tiếp  cận  đối  ngẫu  trong  phương

pháp tuyến tính hóa tương đương áp dụng trong bài toán phân tích đáp ứng của

thiết diện cánh chịu lực khí động.

-  Xây dựng  các  cải  tiến  có  hiệu quả  cho  tiêu chuẩn  tuyến  tính  hóa  tương

đương đối ngẫu, áp dụng cho bài toán flutter phi tuyến của thiết diện cánh.

2

-  Tăng độ chính xác cho nghiệm của bài toán ổn định flutter thiết diện cánh

bằng cách áp dụng tiêu chuẩn đối ngẫu được cải tiến.

-  Thu  được  các  kỹ  thuật  tính  toán  theo  tiêu  chuẩn  đối  ngẫu  cho  bài  toán

flutter của thiết diện cánh.

Thiết diện cánh máy bay theo mô hình hai chiều chịu tác động của lực khí

3. Đối tượng nghiên cứu

động.

4. Nội dung nghiên cứu

4.1. Phương pháp nghiên cứu

- Sử dụng các phương pháp của cơ học để xây dựng mô hình tính toán. Áp

dụng  lý  thuyết  khí  động  học  xây  dựng  phương  trình  dao  động  flutter  của  mô

hình thiết diện cánh máy bay.

- Sử dụng các phương pháp giải tích, đặc biệt phát triển phương pháp tuyến

tính hóa tương đương của cơ học phi tuyến.

- Sử dụng các phương pháp CFD, phương pháp số mô phỏng hệ phi tuyến,

các số liệu thực nghiệm đã có để so sánh, đánh giá kết quả lý thuyết.

4.2. Hướng giải quyết

Trên cơ sở hoàn thiện mô hình cơ học và các kết quả lý thuyết đã có về dao

dựng một kỹ thuật tính toán mới với cách mở rộng tiêu chuẩn đối ngẫu cho bài

toán phân tích đáp ứng phi tuyến của thiết diện cánh dưới tác động của lực khí

động.

động của thiết diện cánh đề tài tập trung phát triển cách tiếp cận đối ngẫu để xây

4.3. Kết quả dự kiến

-  Xây dựng thành công các cải tiến có hiệu quả cho tiêu chuẩn đối ngẫu

cho hệ dao động phi tuyến tuần hoàn và ngẫu nhiên.

-  Áp dụng cho mô hình dao động 2 chiều của thiết diện cánh, xác định các

hiện tượng mất ổn định flutter và vận tốc gió tới hạn.

3

-  Đánh giá sai số của nghiệm và đảm bảo sai số của nghiệm được cải thiện

so với các kết quả đã có trước đây.

5. Cấu trúc của luận án

Luận án gồm phần Mở đầu, 4 Chương và Kết luận.

Mở đầu.

Phần mở đầu trình bày tính cấp thiết, mục đích và nhiệm vụ nghiên cứu, ý

nghĩa khoa học và ý nghĩa thực tiễn của luận án.

Chương 1.

Trong chương này luận án trình bày các kiến thức cơ sở liên quan đến lĩnh

vực khí đàn hồi, sự tương tác giữa ba loại lực: khí động, đàn hồi và quán tính.

Đã tổng quan các nghiên cứu quốc tế và trong nước liên quan đến bài toán phân

tích đáp  ứng của  thiết diện  cánh  chịu lực khí  động.  Các  vấn đề cơ  bản về mô

hình  hóa  thiết  diện  cánh,  các  hiện  tượng  phi  tuyến  và  cách  tiếp  cận  đối  ngẫu

được trình  bày  nhằm  làm  sáng  tỏ  vấn  đề nghiên cứu.  Qua  đó  đã xác định  nội

dung cơ bản cũng như các giới hạn nghiên cứu của luận án.

Chương 2.

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 chuyển động trong dòng khí không nén

được. 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. Sau khi thiết lập phương trình dao động hai

toán số để phân tích các hiện tượng dao động flutter.

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 và tính

Chương 3.

Trình bày tiêu chuẩn đối ngẫu có trọng số cho vấn đề tuyến tính hóa tương

đương hệ dao động phi tuyến, trong đó khi cho tham số trọng số bằng không sẽ

thu được tiêu chuẩn tuyến tính hóa kinh điển. Để giải quyết vấn đề chọn giá trị

tham số trọng số như thế nào sẽ nghiên cứu đề xuất 3 cách lựa chọn tương ứng

với 3 cải tiến. Áp dụng cho các hệ phi tuyến dạng đa thức là dạng hay gặp trong

4

bài toán ổn định flutter 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.

Chương 4.

Tác  giả  sử  dụng  các  phương  pháp  số  cho  phương  trình  vi  phân  chuyển

động của thiết diện cánh. Các kết quả mô phỏng số tính toán đáp ứng phi tuyến

cho  thiết  diện  cánh  máy  bay  được  thực  hiện.  Kết  quả  của  phương  pháp  mô

phỏng số và các kết quả của tác giả khác sẽ dùng để đánh giá, so sánh các kết

quả giải tích thu được theo kỹ thuật tính toán đối ngẫu.

Kết luận chung.

Danh  sách  công  trình  đã  được  công  bố,  đã  được  chấp  nhận  và  sẽ  đăng

thuộc 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

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.

Mechanics, vol. 37, N3, 2015, p.217:230.

5

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, (Accepted 10/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, (Accepted 10/2016).

6

CHƯƠNG 1. TỔNG QUAN VỀ BÀI TOÁN PHÂN TÍCH ĐÁP ỨNG CỦA

THIẾT DIỆN CÁNH

1.1. Khái niệm cơ bản về khí đàn hồi

Khí đàn  hồi  (Aeroelasticity)  là ngành  khoa học  nghiên  cứu  các  hiện  tượng

xảy ra do sự tương tác giữa lực khí động (aerodynamic), lực quán tính (inertia)

và lực đàn hồi (elastic). Lĩnh vực nghiên cứu này được tóm tắt rõ ràng nhất bởi

Lực quán tính

Dao động

Ổn định và  điều khiển

Khí đàn hồi  động

Lực khí động

Lực đàn hồi

Khí đàn hồi  tĩnh

tam giác khí đàn hồi Collar (Collar, 1978) [26], cho trên Hình 1.

Hình 1. Tam giác khí đàn hồi của Collar

Hình 1 thể hiện  mối tương quan của những lĩnh vực như  ổn định và điều

khiển  (stability  and  control),  dao  động  (vibration)  và  khí  đàn  hồi  tĩnh  (static

aeroelasticity) với các tương tác giữa 2 trong 3 loại lực. Sự tương tác của cả 3

loại  lực  này  dẫn  tới  đối  tượng  nghiên  cứu  gọi  là  khí  đàn  hồi  động  (dynamic

aeroelasticity).

Các hiện tượng khí đàn hồi có ảnh hưởng lớn tới việc thiết kế và hiệu quả

hoạt động  của  máy bay.  Lịch  sử  phát  triển của khí  đàn hồi  và  ảnh  hưởng của

môn khoa học này tới việc thiết kế máy bay có thể được tham khảo trong các tài

7

liệu  của  Collar  (1978)  [26],  Garrick  và  Reid  (1981)  [37],  Flomenhoft  (1997)

[34], với những khảo sát về các ứng dụng được cho bởi Friedmann (1999) [35],

và Livne (2003) [62]. Các tài liệu tổng kết khá toàn diện về khí đàn hồi gần đây

bao  gồm  Hodges  và  Pierce  (2002)  [42],  Dowell  vcs  (2015)  [32],  trong  đó  các

cách tiếp cận toán học sâu sắc và các nhiều khía cạnh cơ bản đã được trình bày

chi tiết.

Khí đàn hồi không chỉ là lĩnh vực thuần túy liên quan đến máy bay. Đề tài

này còn rất liên quan tới thiết kế các kết cấu như cầu, xe đua công thức 1, cánh

quạt tua bin gió, cánh quạt động cơ tu bô, máy bay trực thăng, và rất nhiều các

ứng  dụng  khác  …  Trên  thực  tế,  các  nguyên  lý  cơ  bản  cho  các  nghiên  cứu  về

cánh máy bay đều có thể liên quan tới các ứng dụng trên. Các ứng dụng đó đang

ngày  càng  tăng  lên  về  số  lượng  vì  công  nghệ  trong  các  lĩnh  vực  này  đòi  hỏi

những  kết  cấu  nhẹ  hơn  nhưng  làm  việc  trong  điều  kiện  dòng  chất  lưu  khắc

nghiệt  hơn. Các  vấn  đề  này có  thể  được  tham  khảo trong  cuốn  sách  mới nhất

tổng hợp các bài giảng  về khí đàn hồi (Dowell vcs 2015) [31],   và các tài liệu

được trích dẫn trong đó.

Trong luận án này ta tập trung vào hiện tượng khí đàn hồi động (là tâm của

tam giác khí đàn hồi trong Hình 1). Khí đàn hồi động liên quan tới hiệu ứng dao

động của sự tương tác khí đàn hồi, và lĩnh vực chính cần quan tâm là hiệu ứng

phá  hủy  thảm  khốc  của  hiện  tượng  mất  ổn  định  flutter.  Sự  mất  ổn  định  này

không có lợi giữa ba loại lực: khí động, quán tính và đàn hồi, trong đó kết cấu có

thể hấp thụ rất mạnh năng lượng từ dòng khí và bị phá hủy do dao động tăng đột

biến.

thường  liên  quan  tới  hai  hay  nhiều  dạng  dao  động  và  sinh  ra  do  sự  kết  hợp

1.2. Các nghiên cứu đáp ứng của thiết diện cánh

Phân tích đáp ứng của thiết diện cánh là bài toán quan trọng phục vụ quá

trình thiết kế, chế tạo, vận hành và bảo dưỡng máy bay. Để tăng lực nâng, giảm

lực cản, cánh cần được thiết kế có dạng mỏng. Điều này lại dẫn tới độ nhạy cảm  8

với dao động tăng lên, đặc biệt khi chuyển động trong dòng khí với vận tốc lớn.

Lúc này, lực khí động (lực nâng) tăng rất lớn làm biến dạng hình học của cánh,

từ đó lại làm thay đổi đặc trưng của dòng khí, dẫn tới các hiện tượng tương tác

khí đàn hồi. Vấn đề khoa học này đã thu hút nhiều nghiên cứu trong những thập

niên  trở lại đây,  nhất là nghiên  cứu phục vụ nhu cầu chế tạo các loại máy bay

với nhiều tính năng, đảm bảo ổn định khi bay ở các độ cao, vận tốc và điều kiện

bay khác  nhau.  Về mặt mô hình hóa,  cánh có thể được mô  tả từ  đơn giản đến

phức tạp, ví dụ như cho trên các hình 2-5.

Hình 2. Mô hình thiết diện cánh hai chiều, lực khí động quy về lực tập trung

Hình 3. Mô hình cánh theo tấm bị ngàm

9

Trên Hình 2, cánh được mô tả đơn giản là một mặt cắt điển hình, Hodges

và Pierce (2002) [42], Dowell vcs 2015 [31]. Đó là một mô hình vật lý đơn giản,

có phần được sắp đặt trước, rất hữu ích cho việc mô tả các hiện tượng khí đàn

hồi.  Đó  là  sự  đơn  giản  hóa  của  một  hệ  khí  đàn  hồi  bao  gồm  một  tấm  cứng,

phẳng có mặt cắt dạng hình lá (airfoil) được gắn trên hai lò xo đính vào tường

của  đường  hầm  gió.  Trên  mô  hình  đơn  giản  này,  các  nghiên  cứu  lý  thuyết  và

thực nghiệm có thể được thực hiện một cách hiệu quả, đối với cả các hệ tuyến

tính và phi tuyến. Người ta thấy  rằng vận  tốc tới hạn trong hiện  tượng mất ổn

định flutter tính từ mô hình thiết diện cánh hai chiều này có thể xấp xỉ vận tốc

tới hạn của  cánh trong thực tế (Fung 1993) [36]. Nói chung,  vị trí của mặt cắt

điển hình  có thể được chọn ở khoảng 0.7 lần  sải  cánh  tính  từ  gốc.  Mô  hình  2

chiều này cũng rất thích hợp cho việc nghiên cứu thí nghiệm các hiện tượng phi

tuyến và vấn đề điều khiển (Strganac vcs 2000) [78]. Các tham số phi tuyến, các

đáp ứng của mặt cắt có thể được đo trực tiếp. Các tham số của mặt cắt có thể dễ

dàng được điều khiển để khảo sát sự ảnh hưởng. Chính vì những lý do trên, mặc

dù mô hình trên Hình 2 khá đơn giản nhưng vẫn rất hấp dẫn các nhà nghiên cứu

cả về mặt lý thuyết lần thực nghiệm. Mô hình này có thể thực tế hơn nếu xem

xét các điều kiện biên tại đầu cánh. Chẳng hạn như trên Hình 3, cánh được mô tả

bởi tấm được ngàm tại 1 đầu (Wright vcs 2007) [87]. Về cơ bản, mô hình này

cũng  dẫn tới những hiện tượng khí  đàn hồi giống như  mô hình cánh hai chiều

Sự phát triển của máy tính làm tăng khả năng mô hình hóa bằng máy tính

điện tử. Bản thân bài toán khí đàn hồi là sự tương tác qua lại giữa một kết cấu

đàn hồi và môi trường chất lỏng (khí) bao quanh. Do đó , dẫn tới sự phát triển

của  các  phương  pháp  CFD  (computational  fluid  dynamics)  và  CSM

(computational structural model) cũng như sự kết nối giữa chúng (Henshaw vcs

2007) [41].

Hình  4  mô  tả  phương  pháp  CFD  được  sử  dụng  để  tính  toán  các  lực  tác

động lên kết cấu cánh cứng hai chiều.

trên Hình 2, nhưng các công thức thu được sẽ phức tạp hơn.

10

Hình 4. Mô hình cánh phẳng, lực khí động tính bằng CFD

Hình 5. Mô hình kết hợp CFD-CSM

được sử dụng để mô tả kết cấu máy bay đàn hồi còn các phương pháp CFD được

sử dụng để mô tả dòng khí. Sự kết nối của 2 loại phương pháp này là một hướng

nghiên cứu lớn vẫn đang được phát triển mạnh, thường được gọi là sự tương tác

chất lỏng-kết cấu (fluid structure interaction).

Khi  so  sánh  các  mô  hình  ta  thấy  hiển  nhiên  là  mô  hình  phức  tạp  hơn  sẽ

chính xác hơn, tuy nhiên việc khảo sát ảnh hưởng của các tham số cũng như việc

làm nổi bật bản chất vật lý và bản chất toán học của hiện tượng sẽ khó khăn hơn.

Hình 5 thể hiện một mô hình phức tạp hơn trong đó các phương pháp CSM

11

Để  mở đầu cho những nghiên cứu vấn đề tương tác khí động lực theo cách tiếp

cận đối ngẫu, trong luận án này ta giới hạn xem xét theo mô hình trên Hình 2,

trong đó mặt cắt điển hình được đỡ bởi các lò xo phi tuyến và các lực khí động

được tính toán từ lý thuyết cánh mỏng, có thế được quy về lực tập trung tại tâm

khí động (Fung 1993) [36].

Điều khiển máy bay cũng là một vấn đề thu hút được rất nhiều nghiên cứu.

Điều khiển máy bay thông thường được thực hiện bằng cánh nhỏ (aileron) cho

chuyển  động  cuộn  (roll),  bằng  đuôi  lái  cho  chuyển  động  rẽ  (yaw)  hoặc  bằng

bánh lái độ cao (elevator) cho chuyển động lên xuống (pitch). Điều khiển bằng

bánh  lái  nhỏ  cũng  có  thể  được  sử  dụng  để  hạn  chế  hiện  tượng  mất  ổn  định

flutter. Mô hình đơn giản thiết diện cánh điển hình 2 chiều với lực khí động tập

trung  như trên Hình 2 vẫn có thể sử  dụng tiếp tục  đối với bài toán điều khiển

bằng cánh nhỏ, Hình 6.

Hình 6. Mô hình thiết diện cánh 2 chiều có cánh nhỏ điều khiển

Vì mục đích chính của luận án là nghiên cứu cách tiếp cận đối ngẫu cho bài

toán đáp  ứng của thiết  diện  cánh  nên  sẽ không  đi  sâu  vào  các  thuật  toán điều

khiển.  Thay vào  đó,  luận  án  chỉ  xem  xét bài  toán điều  khiển  đơn giản  là điều

khiển PID và khảo sát hiệu quả của cách tiếp cận đối ngẫu trong việc khảo sát

ảnh hưởng của các tham số đến điều khiển PID.

12

1.3. Thiết diện cánh phi tuyến

Phương trình chuyển động ứng với mô hình thường là hệ tự dao động và có

tính chất phi tuyến. Nói chung, khi nghiên cứu hệ phi tuyến này, người ta quan

sát thấy hiện tượng mà ở đó có xuất hiện vòng giới hạn, các hiện tượng rẽ nhánh

Hopf và hiện tượng flutter.  Trong nghiên  cứu  của  mình, Yang  và Zhao (1988)

[89] đã thực hiện cả các tính toán lý thuyết và kiểm chứng thực nghiệm cho mô

hình cánh hai chiều và thu được các kết quả khá phù hợp với nhau. Nghiên cứu

tiếp theo về vòng giới hạn được thực hiện bởi Liu và Zhao (1992) [56]. Các tác

giả  sử  dụng  phương  pháp  cân  bằng  điều  hòa  nhằm  chỉ  ra  những  thông  tin  đủ

chính xác của hiện tượng rẽ nhánh của thiết diện cánh khi có sự thay đổi của tốc

độ dòng khí. Hai ông còn tìm kiếm các nghiệm giải tích dựa vào phương pháp

trung bình và phương pháp khai triển tiệm cận kết hợp với lý thuyết về đa tạp

trung tâm. Tiếp đó một mô hình cánh khí động với tính phi tuyến kết cấu được

nghiên cứu trong bài báo của Yang (1995) [88]. Ông cũng chỉ ra vòng giới hạn

xuất hiện trong mô hình phi tuyến, kết quả này thu được từ việc sử dụng phương

pháp tuyến tính hóa tương đương dựa vào nghiệm tiệm cận bậc hai của phương

pháp  tiệm  cận  của  Krylov-Bogoliubov-Mitropolskii.  Mới  đây,  kết  quả  nghiên

cứu về vòng giới hạn được làm phong phú thêm từ công trình của Shahrzad và

Mhazoon (2002) [76]. Các tác giả đã dự báo rẽ nhánh Hopf trên cơ sở phương

pháp cân bằng điều hòa và phương pháp đa tạp trung tâm. Đóng góp của công

flutter với hai mô hình tuyến tính và phi tuyến có xét đến tính chất khí động ổn

định/không ổn định của hệ.

Nghiên  cứu  đáp ứng  phi  tuyến  của  thiết diện  cánh  khi xảy  ra  hiện  tượng

flutter  cũng  được  nhiều  tác  giả  quan  tâm.  Chẳng  hạn  công  trình  của  Ding  và

Wang  (2006)  [28] đã  nghiên  cứu  hiện  tượng  flutter  khi  vận  tốc  của  thiết  diện

cánh đạt giá trị trên âm. Ổn định của hệ tuyến tính hóa được phân tích trong lân

cận điểm cân bằng. Các tác giả chỉ ra sự mất ổn định flutter là kết quả của hiện

tượng rẽ nhánh  Hopf.  Gần đây  nhất,  Chen  Feixin  và  cộng  sự  phân  tích  flutter

trình này là đã chỉ ra các dao động của thiết diện cánh sau khi xảy ra hiện tượng

13

của thiết diện cánh sử dụng phương pháp tuyến tính hóa tương đương (Chen vcs

2013) [22].  Ngoài ra, các  vấn đề  về  tính toán  số  cho  bài  toán  thiết diện  cánh

được trình bày trong một nghiên cứu của Lee và cộng sự (Lee 1997) [50]. Một

nghiên cứu khá đầy đủ về bài toán thiết diện cánh có thể được tìm thấy trong bài

báo tổng kết (Lee vcs 1999) [51].

Như  vậy  có  thể  thấy  rằng  các  nghiên  cứu  về  thiết  diện  cánh  máy  bay

thường tập trung tính toán đáp ứng nhằm chỉ ra vòng giới hạn trong bài toán phi

tuyến, các hiện tượng rẽ nhánh Hopf và hiện tượng flutter. Phương pháp mà các

tác  giả  sử  dụng  là  phương  pháp  cân  bằng  điều  hòa,  phương  pháp  tiệm  cận,

phương  pháp  tuyến  tính  hóa  tương  đương  và  phương  pháp  đa  tạp  trung  tâm.

Những  đóng  góp  mới  về  các  phương  pháp  này  đã  được  chỉ  ra  trong  từng  bài

toán cụ thể của thiết diện cánh ở trên.

Trên  thực tế,  tính  phi tuyến  của các hệ  khí  đàn  hồi bao  gồm các  tính  phi

tuyến  kết  cấu  và  tính  phi  tuyến  khí  động  lực.  Các  hiệu  ứng  khí  động  lực  phi

tuyến rất khó xác định và đòi hỏi các phương pháp CFD với khối lượng và thời

gian  tính  toán  khá  lớn,  (Djayapertapa  vcs  2001a,  2001b,  Allen  vcs  2005)

[29,30]. Tính phi tuyến kết cấu có thể phát sinh từ sự mòn bản lề của các bề mặt

điều khiển, sự lỏng của các mối liên kết điều khiển cũng như ứng xử phi tuyến

của vật liệu và các nguồn phát sinh khác. Tính phi tuyến kết cấu có thể phân bố

hoặc tập trung tại một vài vị trí. Tính phi tuyến tập trung thường xuất hiện trong

vcs 2003) [51] [32]. Với mục đích nghiên cứu cách tiếp cận đối ngẫu,  luận  án

cũng  sẽ chỉ tập trung vào tính phi tuyến kết cấu tập trung.  Hai dạng điển  hình

của phi tuyến kết cấu có thể được thấy trên Hình 7.

các cơ cấu điều khiển hoặc các bộ phận kết nối của cánh (Lee vcs 1999, Dowell

14

c ụ h p   i ồ h   c ự L

c ụ h p   i ồ h   c ự L

Dịch chuyển

Dịch chuyển

Bậc 3

Khe hở tự do

Hình 7. Phi tuyến bậc ba và phi tuyến khe hở tự do

Sự kết hợp giữa phi tuyến khe hở tự do và phi tuyến bậc 3 được nghiên cứu

bởi Zhao và Hu (2004) [90]. Dạng phi tuyến đa thức hữu tỷ có thể được sử dụng

để mô tả các quan hệ phi tuyến. Mô hình đa thức xấp xỉ có nhiều thuận lợi để

mô tả quan hệ phi tuyến bằng cách điều chỉnh các tham số dựa vào các số liệu

đo. Đa thức xấp xỉ có tính liên tục và khả vi, do vậy phù hợp với nhiều cách tiếp

cận lý thuyết.

Nhiều kỹ thuật được sử dụng trong phân tích hệ phi tuyến có thể được kể

đến bao gồm phương pháp cân bằng điều hòa (Lee vcs 2005, Liu vcs 2005) [52],

phương pháp độ cản tương đương (Chen 2011) [24], phương pháp cân bằng điều

hòa bậc cao, lý thuyết đa tạp trung tâm (Liu vcs 2000) [59], phương pháp ánh xạ

điểm (Liu vcs 2002) [60], phương pháp nhiễu (Chung vcs 2007) [25] và phương

pháp tuyến tính hóa tương đương (Chen vcs 2012, 2013) [22,23]

Các  phương  pháp  được  đề  cập  ở  trên  (trừ  phương  pháp  tuyến  tính  hóa

tự  do  lớn  hơn thì  thông thường  cần  dùng các phương pháp  số,  chẳng  hạn  như

giải phương trình vi phân bằng phương pháp Runge Kutta kinh điển (Tang vcs

2004a) [81]. Ngoại trừ  việc bay thử  nghiệm, các thí nghiệm hầm gió  là những

phương thức có thể xác nhận các kết quả lý thuyết và mô  phỏng số (Tang vcs

2004b, 2000, Strganac 2000) [82] [78]

tương đương) thường chỉ có thể áp dụng cho các hệ có ít bậc tự do. Khi số bậc

15

Hiệu  ứng  chính  của  hiện  tượng  phi  tuyến  là  các  dao  động  vòng  giới  hạn

(limit  cycle  oscillations:  LCOs),  có  thể  được  xem  như  hiện  tượng  dao  động

flutter bị chặn và một ví dụ của loại dao động này cho trên Hình 8.

Hình 8. Dao động vòng giới hạn điển hình

Đôi khi hiện tượng này được nhắc đến với tên gọi là hiện tượng flutter phi

tuyến.  Có  thể  giải  thích  điều  này  trong  trường  hợp  độ  cứng  phi  tuyến  bậc  3

(Hình 7) như sau. Khi vận tốc đạt tới một giá trị nào đó phụ thuộc vào độ cứng

tại vị trí không biến dạng, hiện tượng flutter bắt đầu xảy ra và chuyển động mất

ổn  định  bắt  đầu  hình  thành.  Tuy  nhiên,  khi  biến  dạng  trở  nên  lớn  hơn  thì  độ

cứng cũng trở nên lớn hơn và chuyển động khi đó sẽ bị giới hạn lại chứ không

tiến  tới  vô  hạn.  Trong  một  vài  trường  hợp,  vòng  giới  hạn  được  tạo  thành  bởi

nhiều dao động điều hòa. Rất nhiều các nghiên cứu được thực hiện để tìm kiếm

các phương thức dự đoán vòng giới hạn (LCO) một cách chính xác và hiệu quả.

Hiện tượng LCO cũng chính là đối tượng nghiên cứu chính của luận án.

Tại Việt Nam cũng có một số nghiên cứu về ảnh hưởng của lực khí động

nói  chung cũng  như  hiện  tượng  flutter  trong  nhiều  vấn  đề  kỹ  thuật  khác  nhau

như hàng không, cầu dây văng nhịp lớn, các kết cấu tấm, vỏ,…Đối với vấn đề

khí  động  lực  của  cánh  có  thể  nhắc  tới  một  số  công  trình  sau.  Tác  giả  Lã  Hải

Dũng, (2005) [2] đã khảo ảnh hưởng của một số thông số và nhận thấy nếu tiếp

tục tăng tốc độ bay, dao động của cánh trở thành dao động tự kích gọi là flutter

uốn- xoắn cánh. Việc giải bài toán trên sẽ đưa về giải hệ phương trình vi phân

1.4. Một số nghiên cứu liên quan ở trong nước

16

dao động để tìm ra các thành phần chuyển vị uốn và xoắn cánh theo thời gian ở

các tốc độ bay khác nhau. Trong báo cáo tại Hội nghị Cơ học toàn quốc, tác giả

Hoàng Thị Bích Ngọc vcs,  (2009) [1], đã 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. Các vấn đề dao động của kết cấu cầu dưới tác

dụng của lực khí động cũng được nhiều tác giả Việt Nam quan tâm nghiên cứu.

Trong công trình của Phạm Duy Hòa vcs (2014) [6], đã đặt vấn đề kiểm soát ổn

định flutter của kết cấu cầu hệ treo bằng cách thay đổi hình dáng của mặt cắt cầu

và các khe slot.  Qua khảo sát đã thu được kết quả là vận tốc tới hạn có thể tăng

lên đáng kể bằng cách thay đổi sự tác động của dòng khí lên mặt cắt cầu. Trần

Thế Văn (2012) [8] đã nghiên cứu ổn định của tấm composite lớp chịu tải trọng

khí động bằng phương pháp phần tử hữu hạn có xét đến các phương pháp giảm

dao động bằng thiết bị tiêu tán năng lượng. Trong luận án tiến sĩ của Trần Ngọc

An (2014) [7] đã phát triển phương pháp bước lặp của M. Matsumoto tính vận

tốc gió tới hạn của mặt cắt dầm chủ 3 bậc tự do sang tính toán mô hình mặt cắt

cầu có lắp bộ điều chỉnh rung 4 bậc tự do.  Tác giả đã trình bày nhận dạng tác

dụng của gió và mô hình dao động flutter của dầm chủ trong kết cấu cầu hệ dây,

tính toán ổn định flutter của dầm chủ cầu treo theo mô hình mặt cắt hai bậc tự do

bằng phương pháp bước lặp. Ngoài ra trong luận án cũng tính toán điều khiển

thụ  động dao động  flutter của dầm  chủ  cầu  treo  bằng  phương  pháp  cơ  học  và

bằng phương pháp khí động.

tính  tương  đương  còn  khá  mới  ở  Việt Nam. Điển  hình  là luận  án  của Nguyễn

Ngọc Linh (2015) [3] với tiêu đề “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”.  Đây là  luận án  đầu  tiên  phát  triển

xây dựng các biểu thức cho tiêu chuẩn đối ngẫu và tiêu chuẩn đối ngẫu có trọng

số theo phương pháp trung bình bình phương tối thiểu dựa trên quan điểm đối

ngẫu trong bài toán tuyến tính hóa tương đương. Tiêu chuẩn đối ngẫu có trọng

số được coi như một dạng mở rộng của tiêu chuẩn đối ngẫu và tiêu chuẩn kinh

điển, áp dụng cho quá trình dừng Gauss trung bình không. Dựa trên ý nghĩa của

Đối với các nghiên cứu về cách tiếp cận đối ngẫu cho phương pháp tuyến

17

mức độ phụ thuộc tuyến tính và ảnh hưởng của các quá trình thay thế lượt đi và

lượt về, đã đề xuất phân loại mức độ phụ thuộc tuyến tính của hàm phi tuyến so

với hàm tuyến tính tương đương và phân tích được các đặc điểm và tính chất cơ

bản  của  các  tiêu  chuẩn  này.  Ngoài  ra,  đã  xây  dựng  tham  số  trọng  số  là  hàm

tuyến tính từng đoạn của mức độ phụ thuộc tuyến tính, trong đó dạng giải tích

của trọng số được xác định trên cơ sở đặc điểm của các tiêu chuẩn đối ngẫu, đối

ngẫu có trọng số và bài toán nội suy từ trọng số chính xác của dao động đàn hồi

phi tuyến Lutes Sarkani.

Sự khác biệt cơ bản giữa nội dung luận án của Nguyễn Ngọc Linh với nội

dung nghiên cứu trong luận án này là luận án của Nguyễn Ngọc Linh giới hạn

trong các dao động ngẫu nhiên phi tuyến, trong khi trong luận án này chủ yếu là

dao động tuần hoàn. Ngoài ra trong hai luận án cũng khác nhau ở việc xác định

giá trị trọng số theo các cách tiếp cận khác nhau.

1.5. Cách tiếp cận đối ngẫu

Như đã trình bày ở trên, mặc dù có nhiều phương pháp đã được đề xuất và

nghiên cứu, phương pháp tuyến tính hóa tương đương vẫn là một phương pháp

được ưa thích vì tính đơn giản và hiệu quả trong nghiên cứu các hệ phi tuyến.

Đối với các hệ nhiều bậc tự do thì phương pháp tuyến tính hóa tương đương tỏ

ra  hiệu  quả  hơn  hẳn  các  phương  pháp  phức  tạp  khác.  Tuy  nhiên,  trong  nhiều

trường  hợp,  các  phương  trình  chuyển  động  của  thiết  diện  cánh  có  thể  là  phi

đạt được độ chính xác mong muốn. Mới đây, một cách tiếp cận mới cho bài toán

phi  tuyến  về  dao  động  và  điều  khiển  kết  cấu  được đề  xuất  bởi  Nguyễn  Đông

Anh (2010) [17]. Cách tiếp cận mới được biết với tên gọi cách tiếp cận đối ngẫu

với quan điểm tạo ra một sự hài hòa trong nghiên cứu, cho phép phát hiện bản

chất  của  vấn  đề  một  cách  đầy  đủ  hơn.  Áp  dụng  cho  lĩnh  vực  tuyến  tính  hóa

tương  đương  đã  dẫn  đến  kỹ  thuật  cực  tiểu  bình  phương  đối  ngẫu.  Ban  đầu

phương pháp được đề xuất trong nghiên cứu dao động ngẫu nhiên của các hệ phi

tuyến với kích động ngoài ồn trắng (Anh vcs 2012a) [13]. Kết quả thu được chỉ  18

tuyến mạnh và phương pháp tuyến tính hóa tương đương kinh điển có thể không

ra rằng với các hệ phi tuyến mạnh, phương pháp cho kết quả khá tốt và phù hợp

với  kết  quả  thu  được  từ  nghiệm  chính  xác,  trong  trường  hợp  hệ  phi  tuyến  có

được  nghiệm  chính  xác,  và  kết  quả  thu  được  từ  nghiệm  mô  phỏng  số,  trong

trường hợp không tìm được nghiệm chính xác của hệ phi tuyến. Tính ưu việt của

cách tiếp cận đối ngẫu sau đó còn được phát triển trong nghiên cứu hệ nhiều bậc

tự do chịu kích động ngẫu nhiên (Anh vcs 2012b) [15]. Ý tưởng của kỹ thuật đối

ngẫu  được mở  rộng  sang  nghiên  cứu  điều  khiển  giảm  dao  động  cho  hệ  TMD

(Anh  vcs 2013)  [14].  Các kết  quả thu được về  đáp ứng  của  hệ  TMD  cũng  tốt

hơn  hẳn  so  với  các  kết  quả  đã  có  trước  đây  (Asami  vcs  1999)  [19].  Mới  đây

trong  luận  án  tiến  sĩ,  Nguyễn  Ngọc  Linh  (2015)  [3]  đã  áp  dụng  và  phát  triển

cách  tiếp  cận  đối  ngẫu  cho  hệ  phi  tuyến  một  bậc  tự  do  chịu  kích  động  ngẫu

nhiên.

Như  đã  thấy  cách  tiếp  cận  đối  ngẫu  có  tính  linh  hoạt  và  có  thể  áp  dụng

được cho nhiều lớp hệ phi tuyến khác nhau. Do vậy luận án có mục đích nghiên

cứu phát triển và áp dụng cách tiếp cận này cho bài toán về thiết diện cánh máy

bay. Nghiên cứu  nhằm  tìm  ra những dáng điệu nghiệm  của bài toán với sai số

nghiệm  nhỏ so  với các nghiệm mô  phỏng  số  trong  trường hợp  hệ đang xét  có

tính  phi tuyến,  thậm  chí  là phi  tuyến  mạnh. Hướng  nghiên cứu  này  chưa  từng

được triển khai cho đến thời điểm hiện nay. Việc triển khai nghiên cứu sẽ tạo ra

khả năng thu được các kết quả mới chính xác hơn các kết quả đã biết, bổ xung

không và vũ trụ.

một cách tiếp cận mới trong nghiên cứu đáp ứng phi tuyến của các kết cấu hàng

1.6. Vấn đề nghiên cứu của luận án

Qua các trình bày tổng quan trong chương 1, có thể thấy rằng bài toán phân

tích đáp  ứng của  thiết  diện  cánh  vẫn  là một  bài  toán  thời  sự, các  bài  báo  liên

quan  vẫn  được  xuất  bản  trong  thời  gian  gần  đây.  Việc  phát  triển  các  phương

pháp giải tích cho phép xác định các lời giải gần đúng đáp ứng thiết diện cánh

đang được quan tâm, trong đó có các nghiên cứu mới đây về phương pháp tuyến

tính hóa tương đương. Luận án đề ra vấn đề nghiên cứu là phát triển cách tiếp  19

cận đối ngẫu cho bài toán phân tích đáp ứng phi tuyến của thiết diện cánh chịu

lực khí động. Do luận án sẽ tập trung vào cách tiếp cận đối ngẫu nên mô hình

thiết diện cánh và đáp ứng của nó sẽ chỉ dừng lại ở mức cơ bản, nhưng vẫn thể

hiện các đặc tính chủ yếu của hiện tượng phi tuyến. Cụ thể luận án sẽ giới hạn

nghiên cứu thiết diện cánh ở các điểm sau:

- Thiết diện cánh được xét là một mặt cắt điển hình, hai chiều.

- Lực khí động tác động vào cánh được quy về lực tập trung theo lý thuyết

cánh mỏng tựa dừng.

- Chỉ nghiên cứu hiện tượng đáp ứng phi tuyến điển hình là dao động vòng

giới hạn LCO.

- Bài toán điều khiển bằng cánh nhỏ sẽ được xem xét với luật điều khiển

đơn giản là luật điều khiển PID.

Các vấn đề nghiên cứu của luận án sẽ bao gồm:

- Phát triển, hoàn thiện các kỹ thuật khác nhau của cách tiếp cận đối ngẫu.

- Áp dụng tính toán đáp ứng cho một số hệ phi tuyến cơ bản.

- Áp dụng tính toán hiện tượng flutter cho thiết diện cánh.

- Đánh giá tính chính xác của các lời giải gần đúng thu được theo kỹ thuật

đối  ngẫu  dựa  trên  sự  so  sánh  với  nghiệm  của  các  tác  giả  khác  và  kết quả  mô

phỏng số.

Trong chương này, luận án đã trình bày các kiến thức cơ sở liên quan đến

lĩnh vực khí đàn hồi, sự  tương tác giữa ba loại lực: khí động, đàn hồi và quán

tính. Đã tổng quan các nghiên cứu quốc tế và trong nước liên quan đến bài toán

phân tích đáp ứng của thiết diện cánh chịu lực khí động. Các vấn đề cơ bản về

mô hình hóa thiết diện cánh, các hiện tượng phi tuyến và cách tiếp cận đối ngẫu

được trình  bày  nhằm  làm  sáng  tỏ  vấn  đề nghiên cứu.  Qua  đó  đã  xác định  nội

dung cơ bản cũng như các giới hạn nghiên cứu của luận án.

Kết luận chương 1

20

CHƯƠNG 2. MÔ HÌNH CƠ HỌC CỦA THIẾT DIỆN CÁNH CHUYỂN

ĐỘNG TRONG DÒNG KHÍ

2.1. Lực khí động dừng và tựa dừng

Lực khí động do dòng khí tác động vào cánh có thể được phân loại từ đơn

giản  đến  phức  tạp  bao  gồm:  dừng  (steady),  tựa  dừng  (quasi-steady)  và  không

dừng (unsteady). Trong luận án này ta chỉ xét lực khí động dừng và tựa dừng với

mục đích minh họa các phương pháp tuyến tính hóa.

2.1.1. Lực khí động dừng

Chiều dài dây cung  thiết diện

Trước hết xét một số thuật ngữ cơ bản của cánh như trên hình 9.

α

Góc tới

Hướng gió

Hình 9. Một số thuật ngữ về cánh

Cánh mỏng là một sự lý tưởng hóa của cánh thực với giả thiết độ dày cánh

bằng  0,  tức  là coi  cánh  như  một đoạn thẳng. Do  đó  với cánh  mỏng  ta thường

quan tâm tới 2 đại lượng quan trọng nhất là góc tới hay góc xung kích (angle of

Ta xét một thiết diện cánh 2 chiều như trên Hình 10.

attack) và chiều dài dây cung cánh (chord length).

Hình 10. Dòng dừng đi qua một thiết diện cánh 2 chiều   21

Giả sử  tâm của  hệ  tọa độ được chọn  tại biên  đầu  của  thiết  diện  cánh  với

trục x nằm dọc theo đường dây cung và trục y vuông góc với nó. Nếu một dòng

khí  2  chiều  với tốc  độ U  và  góc  tới  α  (ở  khoảng  cách  xa  với  thiết  diện  cánh)

chảy qua thiết diện cánh, sự nhiễu loạn xuất hiện trong dòng khí bởi sự tồn tại

của thiết diện cánh  làm cho dòng khí có hướng tiếp tuyến với  thiết diện cánh.

Theo lý thuyết khí động lực (Katz vcs 2001), thiết diện cánh mỏng có thể được

thay  thế  bằng  một  sự  phân bố  liên tục  của  các xoáy  vận  tốc  (dải xoáy: vortex

sheet) trên bề mặt cánh. Nếu gọi cường độ của xoáy trên một đơn vị độ dài sải

cánh và trên độ dài dx theo phương dây cung là (x)dx thì lực nâng tác động lên

phần tử dx có thể được xác định theo định lý Joukowsky như sau:

dL

 

  U x dx

(1)

Và tổng lực nâng xác định bởi:

L



  x dx

c    U 0

(2)

Nếu ta xét chiều dương là khi cánh máy bay nâng mũi lên (cùng chiều đông

hồ) thì tổng mô men do lực nâng tạo ra là âm và được xác định bởi:

M

 

 x dx

c   U x   0

(3)

trong đó c là chiều dài dây cung của thiết diện cánh. Với cánh mỏng thì có thể

giả thiết các xoáy nằm trên trục x.

gọi  là  dòng  khí  có  thế  (potential  flow)  (trừ  duy  nhất  chính  vị  trí  điểm  xoáy)

(Katz vcs 2001) [46]. Hình ảnh dòng khí như vậy được cho trên Hình 11.

Một điểm xoáy tạo ra một dòng khí không xoáy (irrotational flow) hay còn

22

Đường đẳng thế

Đường dòng

Hình 11: Đường dòng và các đường đẳng thế của một xoáy 2 chiều tại gốc tọa độ

Vận tốc sẽ có hướng tiếp tuyến với vòng tròn, tỷ lệ nghịch với khoảng cách

đến tâm xoáy và tỷ lệ thuận với cường độ xoáy. Vì ta có một sự phân bố liên tục

của các điểm xoáy nên vận tốc thực tế là tổng hợp của tất cả các vận tốc do các

điểm xoáy gây ra. Vì sự phân bố là liên tục nên sự chồng chất nghiệm này đưa

c

tới tích phân của thành phần vận tốc trên trục x (trục nằm ngang) có dạng:



  v x

d 

         x

 0 2

(4)

Hệ số 2 xuất hiện là do tích phân đường của vận tốc trên một vòng tròn

bằng cường độ xoáy.

a

v U

Độ dốc của dòng chất lỏng trên thiết diện cánh bằng . Đại lượng này

dốc này bằng 0 và điều kiện biên có dạng:

(5)

a

   0

v U

Như  vậy phân bố  của xoáy được xác định từ  (5).  Thêm vào đó điều kiện

Kutta (c)=0 còn phải được thỏa mãn, với chủ ý là dòng khí rời khỏi đuôi cánh

phải trơn. Bây giờ ta xét phép đổi biến sau (Glauert 1959) [38]:

(6)

x



  1 cos

 

c 2

phải bằng độ dốc của bề mặt thiết diện cánh. Nếu ta xét cánh mỏng phẳng thì độ

23

Khi x thay đổi từ 0 đến с dọc theo dây cung thì  thay đổi từ 0 đến . Phân



bố xoáy được tìm dưới dạng (Glauert 1959) [38]:





cot



sin

 n

A n



 2

n

 1

 U A 2  0 

  

(7)







 . Ta có:

0

 x c

  

Dễ thấy dạng nghiệm (7) thỏa mãn điều kiện Kutta

dx



sin

d 

c 2

(8)



cot



sin

A n

   n sin 

d 







 U A  0 

v



  



 1 n cos

cos



  2    

0

   

Thay (7) và (8) vào (4) ta có:





cos

n





cos

n



d







 

A n

  1 



  1 cos







 1 2

     1   

 U A  0 



cos

cos



 1 n    

0





(9)



Ta sử dụng tích phân rất hữu ích sau (Glauert 1959)



,

n



0,1,2...



cos cos

   n sin  sin

n   d      cos 

0

(10)



Áp dụng (10) vào (9) và rút gọn ta sẽ thu được:







cos

n

A n



n

 1

 v U A  0 

  

(11)



a





cos

 n



0

A 0

A n



n

 1

hay

(12)

   0

A 0

Aa , n

Biểu thức (2), (3), (6),  (8), (12) và (7) được dùng để tính lực nâng và mô

men:

2

(13)

L







a

cot

sin

   a

cU



d

 x dx

 2

c     U 0

  U cU 0

Điều kiện biên (5) khi đó dẫn tới:

24

2

2

M

 



cos





sin

2 c U



d

  x dx



 1 cot

c 2

 2

 1    a 4

c  U x   0

  U U  a 0

(14)

Các biểu thức (13) và (14) thể hiện lực và mô men khí động trong mô hình

2

dừng. Ta thấy mô men tại một điểm có khoảng cách x từ biên đầu cánh bằng:

M









  xL M



x

l e .

c 4

  x 

 cU a    

(15)

Biểu thức (15) cho thấy mô men tại vị trí ¼ dây cung cánh tính từ biên đầu

cánh sẽ bằng 0.  Vị trí này được gọi là tâm khí động  (aerodynamic center) của

cánh. Đây là một kết quả quan trọng của lý thuyết cánh mỏng.

2.1.2. Lực khí động tựa dừng

Với dòng không dừng thì các phương trình cơ bản (1), (4) không còn đúng

nữa và bài toán sẽ trở nên phức tạp. Bài toán không dừng có thể được giải quyết

bằng  công  cụ  ánh  xạ  bảo  giác  (conformal  mapping)  trong  lý  thuyết  hàm  biến

phức (mục 6.7, 13.2, 13.3, 13.4 trong Fung 1993) [36]. Tuy nhiên ở đây ta sẽ xét

bài toán đơn giản hóa bằng cách sử dụng giả thiết tựa dừng (quasi-steady) (Fung

1993). Giả thiết này phát biểu rằng các đặc tính khí động của một thiết diện cánh

chuyển động 2 chiều tại bất kỳ thời điểm nào sẽ bằng các đặc tính khí động của

thiết diện đó trong trường hợp nó chuyển động không đổi với các giá trị vận tốc

tức thời tại thời điểm đó. Giả thiết này sẽ cho phép áp dụng các kết quả của mục

2.2.1 một cách trực tiếp. Khi đó các phương trình (1) và (4) vẫn còn đúng.

và nằm  trên trục x. Xét tấm  có 2 bậc tự do: dịch chuyển thẳng đứng h và góc

xoay a quanh trục tại vị trí x0 sau biên đầu cánh (Hình 12).

Xét một tấm phẳng đặt trong dòng chất lỏng với vận tốc tại vô cực bằng U

25

Hình 12: Dòng không dừng đi qua thiết diện cánh

Giá trị của h là dương nếu hướng xuống dưới và giá trị của a là dương nếu

mũi máy bay hướng lên. Xét hệ tọa độ như trên hình 12. Như cách làm ở mục

trước, tấm phẳng sẽ được thay thế bằng một dải xoáy.

 h  





    x a

x 0

Tại vị trí x của thiết diện cánh, thành phần vận tốc thẳng đứng bằng:

trong đó  h  và a  là các đạo hàm theo thời gian của các bậc tự do. Độ dốc tức

thời  của  thiết diện  cánh  là -a. Vì vận  tốc  chất  lỏng phải  hướng tiếp  tuyến  với

thiết  diện  cánh,  thành  phần  thẳng  đứng  của  vận  tốc  chất  lỏng  phải  thỏa  mãn

x



  a

0x

phương trình:

   a



v U

 h U

 U

(16)

tự như mục trước. Xét phép đổi biến (6) và phân bố xoáy dạng (7), ta xác định

được vận tốc có dạng (11). Thay (11) và (6) vào điều kiện biên (16) ta có:



a





cos

 n

 





  1 cos

 

A 0

A n

x 0



 h U

c 2

 a U

n

 1

  

  

dẫn tới

(17)

  a





,



,



0

n



2,3,...





A 0

x 0

A 1

A n

 h U

c 2

 a U

 a c U 2

  

  

Biểu thức (2), (6), (8), (17) và (7) được dùng để tính lực nâng:

Phương trình (16) tương đương với phương trình (5). Cách thực hiện tương

26

L



  x dx

c    U 0











cot



sin

   sin

d

x 0

 h U

c 2

  a  a c U U 2 2

  

  

  U cU 0

 a  

  

   

   



2  a cU







x 0

 h U

c 3 4

 a U

  

  

  

  

(18)

So sánh (13) và (18) cho thấy lực nâng của mô hình tựa dừng được thay đổi

a







x 0

 h U

c 3 4

 a U

  

  

từ mô hình dừng bằng cách thay a bởi đại lượng . Như vậy

mô hình tựa dừng làm xuất hiện thêm hai thành phần cản khí động.

Mô men đối với trục đi qua vị trí x0 được xác định theo công thức tương tự

a







x 0

 h U

c 3 4

 a U

  

  

2

với (15) nhưng thay góc a bởi đại lượng :

M















x 0

x 0

x 0

c 4

 h U

c 3 4

 a U

  

  a cU  

  

  

  

  

(19)

Các trình bày ở trên dựa trên mô hình cánh mỏng. Trên thực tế, các giá trị

thí  nghiệm  của  lực  nâng  không  thực  sự  trùng  với  lý  thuyết  do  giả  thiết  cánh

mỏng. Vì thế, để tiện lợi cho các phân tích có số liệu lấy từ thí nghiệm, người ta

2

2





viết lực nâng (18) và mô men (19) dưới dạng:

L



C



a







L

c a l

x 0

cU 2

cU 2

 h U

c 3 4

 a U

  

  

  

  

2

2





(21)

M



C









M

c m

a

x 0

2 c U 4

2 c U 4

 h U

c 3 4

 a U

  

  

 a  

  

trong đó CL và CM là hai đại lượng phi thứ nguyên được gọi là hệ số nâng và hệ

số mô men, cla và cma là hai đại lượng phi thứ nguyên là hệ số nâng và hệ số mô

men trên một đơn vị góc. Theo lý thuyết (công thức  (18), (19)) thì:

(22)



 2 ,





c a l

c m

a

c a l

02 x c

1 2

  

  

(20)

27

Thực tế các mô hình thí nghiệm cho thấy các hệ số này chỉ gần với các hệ

số tính toán trên lý thuyết (Wright vcs 2007) [87].

Trên thực tế, mô hình tựa dừng có thể sử dụng cho các bài toán có tần số

dao động không quá lớn. Trong các bài toán liên quan tới tần số dao động cao

thì mô hình không dừng tổng quát cần phải được sử dụng (Fung 1993) [36].

2.2. Phương trình chuyển động của thiết diện cánh

Xét một dải có độ dày đơn vị của thiết diện cánh phẳng 2 chiều có 2 bậc tự

do: bậc tự do uốn h (chiều dương hướng xuống tính từ trục đàn hồi) và bậc tự do

c = 2b

b

Vị trí đàn hồi

h

a bh

Vị trí dịch chuyển

α

Trục đàn hồi

xoắn a (chiều dương là đầu cánh nâng lên) quanh trục đàn hồi (Hình 13)

Hình 13. Mô hình cánh 2 bậc tự do

được với tốc độ U. Các phương trình chuyển động của thiết diện cánh được xây

dựng từ việc xét cân bằng của  các lực quán tính, lực  đàn hồi và lực  khí động.

Với một phân tố khối lượng dm ở vị trí ở khoảng cách r (chiều dương hướng về

cạnh đuôi) tính từ trục đàn hồi, Lực đàn hồi bằng:



Giả  sử  thiết  diện  cánh  được  đặt  trong dòng  khí  của  chất  lỏng  không  nén

      dm h ra

 mh S 

 





Tổng lực quán tính trên một đơn vị sải cánh bằng:      a

  dm h r

  a





28



 dm

tổng  khối  lượng  của  cánh  trên  một  đơn  vị  sải  cánh  và trong  đó  m

S



rdm





mô men tĩnh của cánh quanh trục đàn hồi. Tích phân được lấy trên

toàn bộ dây cung cánh.

Lực quán tính gây ra một mô men trên một đơn vị sải cánh quanh trục đàn





dm

 



 a I a

hồi, được xác định bởi:

  r h r

  a



  Sh



2 r dm



a  I



trong đó mô men quán tính khối lượng quanh trục đàn hồi.

Giả sử các chuyển dịch uốn và xoắn được đỡ bởi cặp lò xo tại trục đàn hồi

với hệ số lò xo kh và ka, tương ứng. Lực đàn hồi tương ứng với chuyển dịch h

bằng –hkh, theo hướng ngược với h. Lực ngược hướng với a là -aka.

Phương trình chuyển động có thể được viết theo điều kiện là tổng của lực và mô

men quán tính và lực và mô men đàn hồi phải cân bằng với lực khí động bên

 mh S 



 

L

k h h

ngoài (là lực nâng và mô men uốn được xác định từ (20), (21)). Ta có:





M

 a I a

 a  Sh k  a a

   

(23)

Ở đây  ta  chú  ý lực  nâng  được định  nghĩa  là chiều dương  hướng  lên theo

quy ước dấu thông thường trong khí động lực. Trong khi đó chuyển dịch uốn h

lại  được  định  nghĩa  chiều  dương  hướng  xuống.  Do  vậy  lực  bên  vế  phải  của

chú ý rằng khoảng cách từ trục đàn hồi đến biên đầu cánh là

(với



b



x 0

  1h a

b là nửa chiều dài dây cung cánh), ta biến đổi (23) về dạng:

 mh S 

 a



 









k h h

2 bU c a l

a h

 h U

1 2

 a b U

  

  

(24)

2









a

 a I a

 Sh k  a   a

2 b U c m

a

h

 h U

1 2

 a b U

  

  

 a    a  

     

      

phương trình thứ nhất của (23) có dấu trừ. Sử dụng các biểu thức (20), (21) với

29

2.3. Hiện tượng flutter

Một kết cấu có mặt cắt không tròn sẽ chịu lực dòng khí thay đổi theo góc

mà mặt cắt tạo với dòng khí. Khi kết cấu dao động, góc này cũng dao động và

lực  chất  lỏng  cũng dao  động.  Nếu  lực  chất  lỏng  có  xu  hướng  khuếch  đại dao

động thì kết cấu trở nên mất ổn định khí động và biên độ dao động có thể trở lên

rất lớn. Hiện tượng đó gọi là galloping hoặc flutter.

Khái  niệm  galloping  chỉ hiện  tượng mất ổn  định  khí động  xảy  ra  với  các

mặt cắt tù (bluff), tức là các mặt cắt có tách dòng phía sau. Ngược lại, khái niệm

flutter dành cho các mặt cắt khí động, không có tách dòng phía sau.

2.3.1. Hiện tượng mất ổn định 1 bậc tự do

Hiện tượng mất ổn định 1 bậc tự do  một cách tổng quát có thể được giải

thích về mặt nguyên lý như trên hình 14. Lực tác dụng vào mặt chắn phụ thuộc

 arctan y U



vào góc tới a. Góc này bằng trong đó y là chuyển dịch lên xuống

và U là vận tốc dòng khí.  Vậy lực chất lỏng phụ thuộc vào vận tốc dao động.

Trong trường hợp lực này cùng chiều với vận tốc dao động thì sẽ gây ra cản âm.

Với vận tốc đủ lớn thì cản âm này sẽ thắng cản dương của bản thân kết cấu và

gây ra mất ổn định.

Hình 14: Mô hình lực chất lỏng tác động vào hệ 1 bậc tự do

30

Cụ  thể  hơn,  ta  xét  phương  trình  (24).  Nếu  mặt  cắt  chỉ  có  dao  động  lên



xuống, không có dao động xoắn (a=0) thì phương trình dao động (24) còn lại:

 mh k h

 

h

   bUc ha

l

(25)

Trong thực tế, hệ số cla theo lý thuyết bằng 2 (công thức (22)) với các góc

nhỏ. Khi đó hệ (25) có cản dương và luôn ổn định. Tuy nhiên khi góc tới  lớn

hơn một góc giới hạn thì thực nghiệm đã chỉ ra rằng cánh máy bay xảy ra hiện

tượng “chết đứng” (stall). Sau khi đó tăng góc nâng lại làm giảm lực nâng. Đồ

)M C

) D     C

) L     C

n ả c   n ê r t   g n â n   ố s   ỉ

T

(   n ả c   c ự l   ố s   ệ H

(   n e m   ô m   ố s   ệ H

(   g n â n   c ự l   ố s   ệ H

Góc tới (độ)

Góc tới (độ)

thị điển hình của lực nâng theo góc xung kích được cho trên hình 15.

Hình 15. Các hệ số khí động theo góc xung kích

Khi đó hệ (25) có thể có cản dương và xảy ra mất ổn định.

Nếu  mặt  cắt  chỉ  có  dao  động  xoắn,  không  có dao  động  lên  xuống  (h=0),

phương trình dao động (24) còn lại:

2

(26)







 a a  I k a a

2 b U c m

a

a h

 a b U

1   2 

  

 a  

  

Ta thấy rõ ràng từ (26) rằng nếu vận tốc đủ lớn thì

31

2



2 b U c m

 ka

a

(27)

hệ trở thành hệ có độ cứng âm và mất ổn định. Đây còn gọi là mất ổn định tĩnh.

2.3.2. Hiện tượng mất ổn định 2 bậc tự do

Hiện tượng mất ổn định 2 bậc tự do có thể giải thích về mặt nguyên lý như

sau. Trong trường hợp này lực khí động của chất lỏng phụ thuộc vào góc xoắn.

Lực này có thành phần nâng tác động lên chuyển động thẳng đứng của mặt cắt.

Nếu lực nâng có xu hướng cùng chiều với chuyển động thẳng đứng thì dao động

sẽ bị khuếch đại. Ở đây, điều kiện để xảy ra mất ổn định là phải có sự tương tác

giữa 2 dạng dao động: lên xuống và xoắn.  Lực nâng  phụ thuộc vào góc xoắn,

góc xoắn lại bị tác động bởi chuyển động lên xuống. Hình 16 mô tả hiện tượng

mất ổn định kết hợp này.

Hình 16: Mô tả hiện tượng mất ổn định lên xuống –xoắn

Với mô hình dừng mô tả bởi phương trình (24) bỏ qua cản, người ta có thể

đầu xảy ra hiện tượng mất ổn định flutter.

Cụ thể hơn, xét phương trình (24) và bỏ qua các đại lượng liên quan tới vận

tốc, ta có:

 mh S 





 

k h h

2 a bU c a l

(28)

2



 a I a

 a  Sh k  a   a

2 b U c m

a a

   

Đây là một hệ tuyến tính. Xét nghiệm của hệ có dạng có dạng

pt

(29)

 h He

,pt

 Aea

tìm được lời giải giải tích của tốc độ tới hạn, là tốc độ dòng chất lỏng ở đó bắt

32

Trong đó H và A là các biên độ dao động phức. Khi đó, thay (29) vào (28)

S

2

p





0

m  S

k h 0

 

 I a

2 bU c a l 2 2  b U c m

a

k a

  

  

  

  

ta thấy p là nghiệm của bài toán giá trị riêng:

4

2

Hay có dạng:





   0

a p 4

a p 2

a 0

(30)

2





k



,

 a mI 4

a

2 S a ; 0

h

k a

2  c b U m

a

 

với

2

  2







2  c b U m

a

k I h

a

 c b U S a l

 a m k   a 2

 

(31)



4

2

  a 2

a a 4 0

Nghiệm của phương trình có dạng:

p



2 a 2 a 2

4

(32)

Ta  thấy  a4>0  trong  (31)  với mọi phân bố khối  lượng.  Giá trị  của  a0  cũng

dương nếu vận tốc nhỏ hơn vận tốc mất ổn định tĩnh (27). Rõ ràng đây là trường

hợp cần quan tâm vì nếu tốc độ tới hạn flutter 2 bậc tự do nhỏ hơn tốc độ mất ổn

định tĩnh thì ta không cần quan tâm đến nó. Tóm lại là a4 và a0 đều dương.   Ban đầu vận tốc nhỏ thì giá trị của a2 trong (31) là dương. Giá trị của p2 trong



4

   0

2 a 2

a a 4 0

(32) là thực và âm nếu:

nếu xảy ra điều kiện:



4

   0

2 a 2

a a 4 0

thì  giá  trị  của  p2  là phức. Khi đó p  không  thể  là tất  cả thuần  ảo.  Theo định  lý

Viet áp dụng  vào (30) thì tổng  các phần  thực  của  tất  cả của p bằng 0. Khi đó

chắc chắn phải có ít nhất một giá trị của p có phần thực dương. Hệ khi đó mất ổn

định. Lý luận đó dẫn tới điều kiện:

(33)



4

   0

2 a 2

a a 4 0

Khi đó p sẽ là các giá trị thuần ảo, dao động là ổn định cân bằng. Ngược lại

33

là điều kiện biên giữa miền ổn định cân bằng và miền mất ổn định. Phương trình

2

2

2







0

b 4

b 2

b 0

 U 2

 U 2

  

  

(33) là 1 phương trình bậc 4 trùng phương với U:

2

2

2



,







4

mI







b 4

a

Sc a l

b 0

mk a

k I h

a

a



trong đó:

2



mI



2





mk

 







b 2

a

a

h

mbc m

a

Sc a l

a

k I h

a

 b mbc 2 m 

  S c bk m

 

  S k k , a h  b 4

(34)



  b 2

b b 4 0 4

Giải phương trình trùng phương ta có:

 U



2 f

1 2

2 b 2 b 2 4

(35)

Như  vậy  (35)  cho  ta  biểu  thức  giải  tích  của  vận  tốc  tới  hạn  xảy  ra  hiện

tượng mất ổn định flutter trong hệ tựa tĩnh.

Đặc biệt, trong tài liệu (Pines 1958, Dowell vcs 2015) [67] [31],  người ta

đã chứng minh được rằng nếu S≤0, tức là trọng tâm cánh nằm ở trước tâm đàn

hồi thì các nghiệm của (35) đều âm và không xảy ra hiện tượng flutter.

2.4. Tính toán vận tốc flutter trong hệ tuyến tính

Mục này trình bày quy trình chung nhất để xác định vận tốc tới hạn xảy ra

hiện tượng flutter trong hệ tuyến tính tổng quát. Sau đó một tiểu mục của mục

trạng thái, sẽ được sử dụng nhiều trong chương 4.

này được dành để đưa ra công thức cụ thể tính vận tốc tới hạn flutter cho hệ 5

2.4.1. Hệ tự dao động tổng quát

Ta sẽ xét một hệ tự dao động tổng quát được mô tả bởi phương trình trạng

thái:



(36)

z

UMz A



thuộc  vào vận  tốc dòng  khí,  M  là ma  trận  khối  lượng  mở rộng  cỡ  nn .  Hiện

tượng mất ổn định flutter  được  hình  thành  như  sau.  Khi  vận  tốc dòng  khí  còn  34

Trong đó z là vectơ gồm n trạng thái, A là ma trận hệ thống cỡ nn và phụ

nhỏ hơn một giá trị tới hạn thì đáp ứng của (36) là hội tụ về 0 với mọi điều kiện

đầu.  Lúc đó  các giá trị riêng của hệ tuyến tính (36) đều có phần thực  âm. Khi

vận tốc đạt tới giá trị tới hạn, xuất hiện giá trị riêng có phần thực bằng 0, hệ có

thể xuất hiện dao động với biên độ hữu hạn hoặc đáp ứng tiến tới vô cùng.

Ký hiệu giá trị riêng trong trường hợp tới hạn là i trong đó i là số ảo đơn vị còn

 là một số thực dương thể hiện tần số của dao động với biên độ hữu hạn (nếu

nó tồn tại). Vậy i là nghiệm của phương trình đặc trưng:

  i

(37)

M A



0

 U



n

n



Phương trình (37) là một phương trình đa thức bậc n với i dạng:

 i

  ...





   0

 



 1

 a i n

 a

n

1

a i 1

a 0

(38)

Trong đó ai (i=1,,n) là các hệ số thực và phụ thuộc vào tốc độ dòng khí U.



2



 

... 0

a n



n  2

 1

 3



a

q



a

 

... 0

  n    a p n   

n   1

n

n  3

n



Tách phần thực và phần ảo của (38) ta có 2 phương trình đa thức:

Điều  kiện  để  2  đa  thức  này  có  nghiệm  chung  là  định  thức  của  ma  trận

Sylvester (hay kết thức) (xem phụ lục) tương ứng với chúng phải bằng 0. Như

vậy ta có phương trình:

S

R

   0

p q ,

p q ,

(39)

trong đó Sp,q là ma trận Sylvester của hai đa thức, Rp,q là kết thức của 2 đa thức

phương trình để xác định vận tốc tới hạn flutter.

Ngoài ra, nghiệm thu được chỉ có thể tồn tại trong thực tế nếu nó ổn định

với các nhiễu. Việc phân tích ổn định có thể dựa theo tiêu chuẩn sau (Wei vcs

2014) [86]. Giả sử ta đã thu được một dao động vòng tới hạn (LCO) với biên độ

dao động là A, phương trình (38) có nghiệm với phần thực bằng 0. Sự thay đổi

nhỏ của biên độ dao động sẽ dẫn tới sự thay đổi của nghiệm thuần ảo này. Ký

hiệu  sự  thay  đổi  của  phần  thực  của  nghiệm  là  .  Nếu  một  nhiễu  dương  của

(xem  phụ  lục).  Phương  trình  (39)  là  một  phương  trình  với  ẩn  là  U,  chính  là

35

biên độ (ΔA>0) dẫn tới giá trị âm của phần thực (<0) của nghiệm được xét thì

dao động LCO là ổn định bởi vì năng lượng sẽ được tiêu tán cho đến khi biên độ

giảm  về  giá  trị  khi  chưa  bị  nhiễu.  Tương  tự  như  vậy,  một  nhiễu  biên  độ  âm

(ΔA<0) dẫn tới phần thực dương (>0) của nghiệm  cũng tạo ra một LCO ổn

định bởi vì biên độ sẽ được tăng lên cho đến khi đạt được biên độ LCO ổn định

lần nữa. Vậy tóm lại, điều kiện

Δσ.ΔA < 0 (40)

xác định một dao động LCO ổn định.

2.4.2. Thiết diện cánh 2 chiều có điều khiển PID

Trong  chương  4  ta  sẽ  tập  trung  xem  xét  phương  trình  chuyển  động  của

thiết diện cánh 2 chiều có cánh nhỏ được điều khiển bằng thuật toán điều khiển

PID.  Phương trình của hệ này có thể được viết chung dưới dạng phương trình

trạng thái 5 chiều như sau:



1 0 0 0 1 0 0 0 1 0 0 0

0 0 1 k

0 0 0 k

1 0 0 k

0 1 0 k

z 1 z 2 z 3 z

1

2

0 0 0 k 1

2

3

4

5

0 0 0

k

k

2

3

6

7

k 8

k 9

k 10

4 z 5

       

 z 0 0   1    z 0 0   2    0 z 0 3    m m z   4    m m z   5

       

       

               

       

(41)

trong  mi  (i=1,..3)  và  ki  (i=1,..10)  là  các  hệ  số  sẽ  được  xem  xét  chi  tiết  trong

chương 4. Khai triển đa thức đặc trưng, phương trình (38) trong trường hợp cụ

5

4









 

  i

  



2 2

m k 2 9

3

k m m k m k 4

2 5

3

1 10

3











 

 m m m i 1 m k m k m k 1 7

k k 4 10

2 2

2 6

k k 5 9

2











k k 7 4

k k 1 10

k k 2 9











  i 0





  i

       

k k m k m k 6 5 1 8 2 3 

k k 6 2

k k 3 9

k k 8 4

k k 1 7

k k 1 8

 k m i 1 3  k k 6 3

Tách phần thực và phần ảo ta có 2 phương trình đa thức:















 2 

 m m m 3

1

2 2

k k 4 10

m k m k m k 1 7

2 2

2 6

k k 5 9

k m 1 3

 4 

(42)

 









0

k k 6 2

k k 3 9

k k 8 4

k k 1 7

thể này có dạng:

36







 4 

k m m k m k 4

2 5

3

1 10

 

















0

 m k 2 9 

 2 

k k 7 4

k k 1 10

k k m k m k 6 5 1 8 2 3

k k 2 9

k k 1 8

k k 6 3

(43)

giữa 2 phương trình (42) và (43) dẫn tới phương trình để tính tần

Khử 4









k k 4 10

2 2

k k 5 9

k m 1 3





 



1

1 10







2 5 

k k 7 4

k k 2 9

m k 2 9 k k 6 2

k k 1 7





 

 

 

 

 m k m k m k 2 6 1 7 2  m m m 2 3 k k m k m k 6 5 1 8 2 3 k m m k m k 4

k k 1 10 m k 2 9

2 5

3

1 10

 k k k k 6 3 1 8 k m m k m k 4 3  k k k k 3 9 8 4 2  m m m 2 3

1

   2    

     

     

     

số:

(44)

Các  đa thức  (42) và (43) có thể xem  như  các  đa thức bậc 2 với 2.  Điều

2







k k 1 8

k k 6 3

1

2 2

















 k m m k m k 4

2 5

3

1 10

m k 2 9

k k 6 2

k k 3 9

k k 8 4

k k 1 7

   

   



















k k 1 7

k k 3 9

k k 8 4

k k 7 4

k k 1 10

k k 2 9

kiện trùng nghiệm là kết thức của 2 đa thức bằng 0 (xem phụ lục):





  m m m 3     















k k 6 2 

k k m k m k 6 5 1 8 2 3 



k k 4 10

2 2

2 6

k k 5 9

k m 1 3

k k 6 3

k k 1 8

   















k k 1 10

2 2

3

k k m k m k 6 5 1 8 2

3

k k 2 9

















     m m m k k 1 7 4   

 



 m k m k m k 1 7  k m m k m k 4

m k 2 9

2 5

1 10

3

k k 4 10

m k m k m k 1 7

2 2

2 6

k k 5 9

k m 1 3

   

   

(45)

Với  các  giá  trị  cụ  thể  của  mi  (i=1,..3),  ki  (i=1,..10)  trong  các  ví  dụ  của

chương 4, phương trình (45) sẽ là phương trình để xác định vận tốc tới hạn còn

tại).

(44) là phương trình để xác định tần số dao động biên độ hữu hạn (nếu  nó tồn

2.5. Tính toán thiết diện cánh bằng phương pháp CFD

Trong  các  mục  trên,  ta  đã  trình  bày  một số  vấn  đề  lý  thuyết  về  mô  hình

dừng và tựa dừng của dòng khí qua thiết diện cánh. Trong mục này ta trình bày

một số nghiên cứu bằng phương pháp CFD để kiểm tra một số kết quả lý thuyết

về mô hình dừng cũng như tối ưu hóa thiết diện cánh, cải thiện hiệu  suất cánh

máy bay, xác định được góc tới giới hạn. Dựa trên lý thuyết khí động và khí đàn

37

hồi,  các  phân  tích  CFD  thu  được  bằng  cách  sử  dụng  gói  phần  mềm  ANSYS

Fluent. Các mô phỏng đã được tiến hành trên các mô hình tính toán của cánh sử

dụng cho máy bay không người lái.

2.5.1. Mô phỏng khí động lực trên mô hình cánh máy bay

Nhiều  công  cụ mô phỏng  dựa  trên phân  tích  CFD  và  phương  pháp  số đã

được phát triển và có thể chứng minh là cực kỳ hữu ích cho các nghiên cứu về

khí  động  trên  máy bay. Nó  có  tầm quan  trọng  đặc biệt  trong thiết  kế các máy

bay không người lái (UAV) thử nghiệm ở các nước đang phát triển, trong đó có

Việt Nam. Các UAV có thể được sử dụng cho các ứng dụng cả quân sự và dân

sự khác nhau như giám sát ven biển, quan sát thời tiết, theo dõi cháy rừng, thu

thập dữ liệu khoa học, vv.

Lực  khí  động  tác  động  lên  đối  tượng  bay  bao  gồm  lực  nâng  theo  hướng

vuông góc với dịch chuyển bay và lực cản theo hướng  ngược với dịch chuyển

bay. Một loạt các vấn đề về khí động có liên quan đến việc xác định sự tương

tác  giữa  không  khí và  một  vật thể  rắn  di chuyển  trong nó,  như  một  cánh máy

bay.  Các  nhà  nghiên  cứu khí  động  đã quan  tâm  đến  các  tối  ưu  hình  dạng  của

thiết  diện  cánh  máy  bay, để  cung  cấp  lực  nâng  lớn  nhất  và  lực  cản  thấp  nhất

trong quá trình cất cánh, cũng như trong khi bay.

Một  thiết  diện  cánh  được định  nghĩa  như  mặt  cắt của một  cánh  máy bay

(James 1997) [44]. Thiết kế cơ bản của thiết diện cánh máy bay được thể hiện

Biên dẫn

Độ dày

Biên theo

Dây cung

Chiều dài dây cung

trong hình 17.

Hình 17. Phác thảo thiết kế thiết diện cánh máy bay.

Các  nghiên cứu  có hệ thống đầu tiên về hình dạng thiết diện cánh và các

hiệu năng của chúng được xây dựng bởi Ủy ban tư vấn Quốc gia về hàng không,

38

(NACA,  tiền  thân  của  NASA  ngày  nay).  Các  dây cung, độ  cong  và  độ dày  là

những tính năng quan trọng nhất của hình học cánh. Các đặc tính hiệu suất của

các thiết diện cánh thường được bao gồm lực nâng, lực cản, áp suất động lực và

mômen đối  với tâm  khí  động.  Các  giá  trị của  những  đặc  tính  này  đã được đo

bằng các thử nghiệm đường hầm gió và hoặc cũng có thể được xác định thông

qua phân tích lý thuyết toán học, hoặc bằng cách sử dụng mô hình tính toán với

các công cụ mô phỏng CFD.  Trong mục này, phân tích  CFD trên các mô hình

cánh máy bay sử dụng thiết diện NACA 2412 được thực hiện bằng phần mềm

ANSYS. Mục tiêu của nghiên cứu này là thử nghiệm các mô hình cánh chế tạo,

mà có thể được sử dụng để thiết kế các UAV.

Như  đã  trình  bày  trong  các  nghiên  cứu  lý  thuyết,  lực  nâng  theo mô  hình

2

dừng được xác định bởi:

L



C SU

1 2 L

(46)

a 2

cl  là diện tích mặt cánh từ góc chiếu

LC

trong đó là hệ số lực nâng,  S

theo phương lực nâng.

Lực khí động khác vuông góc với lực nâng, được gọi là lực cản. Lực này

chống lại các chuyển động tương đối của cánh máy bay và có hướng song song

với luồng không khí (Prabhakar 2013) [69] xuất hiện do ma sát giữa vỏ với các

phần  tử  không  khí  trên  bề  mặt  của  cánh.  Trong các tính  toán lý thuyết,  ta  xét

cản được tính bởi:

2

(47)

D



C AU

1 2 D

với D là lực cản, A là diện tích hình chiếu tính lực, và CD là hệ số lực cản.

Hệ số lực nâng và hệ số lực cản là hai hệ số không thứ nguyên của các lực

khí động. Các hệ số này có thể xấp xỉ bằng cách sử dụng lý thuyết khí động như

trên, tính toán bằng mô phỏng số hoặc đo đặc trong các thử nghiệm đường hầm

gió của một cấu hình cánh máy bay hoàn chỉnh.

chất lỏng không nhớt nên lực cản tính ra bằng 0. Tuy nhiên, trong thực tế, lực

39

Phần mềm ANSYS cung cấp bộ giải pháp kỹ thuật mô phỏng trong các vấn

đề  kỹ  thuật  mà  một quá  trình  thiết kế  đòi  hỏi.  Các  công  ty  trong  các  lĩnh  vực

công nghiệp thường sử dụng phần mềm này. ANSYS sử dụng mô hình hóa phần

tử hữu hạn (FEM) và nhiều thuật toán lập trình khác cho các mô phỏng và tối ưu

hóa  các  bài  toán  thiết  kế  khác  nhau.  ANSYS  có  nhiều  phần  phụ,  trong  đó

ANSYS  Fluent  và  Structural  đã  được  chọn  để  chạy  các  mô  phỏng  trong  mục

này. CFD được áp dụng để phân tích bài toán về cơ học chất lỏng và động lực.

Các  khả  năng  mô  hình  vật  lý  và  kết  quả  chính  xác,  nhanh  chóng,  cho  thấy

ANSYS là một trong những phần mềm toàn diện nhất cho các mô hình CFD có

sẵn trong thế giới ngày nay.

Một sơ đồ mô hình hình học của cánh máy bay và thiết diện cánh được thể

hiện trong hình 18.

Hình 18. Mô hình cánh máy bay tạo ra trong ANSYS.

Có một số phương pháp được sử dụng  để mô tả hình dạng của cánh máy

NACA  2412  (Ira  1951)  được  sử  dụng  để  thiết  kế  cánh  (Brat  1971)  [20].  Các

thông số được lựa chọn như chiều dài dây cung c=0.3m, chiều dài cánh l=1.6m.

Các  kích  thước này thường được sử dụng trong chế tạo một số mẫu cánh  thực

nghiệm,  có  trong  các  dữ  liệu mở  của  một  số  mẫu  UAV  đã  giới  thiệu  tại  Việt

Nam.

Về vật liệu, luyện kim đã đóng một vai trò quan trọng trong sự phát triển

của ngành hàng không. Cho đến gần đây, một số vật liệu mới đã được áp dụng

trong  xây  dựng  máy  bay,  chẳng  hạn  như  hợp  kim  titan,  hoặc  composite.  Tuy  40

bay, như NACA bốn chữ số, năm chữ số, vv. Trong nghiên cứu này, thiết diện

nhiên, các siêu vật liệu vẫn còn khá đắt đối với sản xuất máy bay, và chưa có các

nghiên  cứu  đầy  đủ.  Với  những  ưu  điểm  về  trọng  lượng  và  chi  phí,  hợp  kim

nhôm vẫn được sử dụng rất rộng rãi.

Hợp kim nhôm 7075 T6 (Ma vcs 2014) [63] cũng là một vật liệu hấp dẫn,

vì vậy nó được sử dụng trong các thiết kế của nghiên cứu này. Các đặc tính vật

liệu liên quan được đưa ra trong bảng 1.

Bảng 1: Các đặc điểm của Hợp kim nhôm 7075-T6

TT Tính chất vật liệu Giá trị Đơn vị

1  Khối lượng riêng  2810 kg/m3

2  Mô đun Young 71.7 GPa

3  Mô đun trượt 26.9 GPa

4  Hệ số Poison 0.33 -

Các  dòng  khí được  xem là  luồng  không khí.  Các đặc  tính dòng  khí được

chọn tương tự như được sử dụng trong các thí nghiệm, chẳng hạn như mật độ là  1.225kg/m3, và độ nhớt dòng khí là 1.7894e-5. Tất cả các thông số của vật liệu

trên được áp dụng để thiết lập cho mô phỏng.

Để phân tích  dòng khí,  các miền dòng khí được chia  thành các miền phụ

nhỏ hơn, được gọi là chia lưới. Mục đích sử dụng của lưới là để tách và tính toán

mô hình không gian. Các lưới sử dụng được hiển thị trong hình 19.

các tính chất của dòng khí không khí. ANSYS Fluent sử dụng các mắt lưới để

Hình 19. Chia lưới trong ANSYS

41

ANSYS Fluent giải quyết các phương trình Navier-Stokes tại mỗi nút của

lưới. Hơn nữa, một phương pháp lặp đi lặp lại được sử dụng bởi ANSYS Fluent

để hội tụ nghiệm trong phân tích này.

Trước  khi  có  thể  chạy  mô  phỏng,  các  thông  số  cụ  thể  và  điều  kiện  biên

được thiết lập. Đầu tiên một số thiết lập chung cần phải được cài đặt. Tiếp theo

các mô hình dòng khí được chọn. Các bước tiếp theo, điều kiện biên được thiết

lập cho các khu vực khác nhau của mỗi lưới, chẳng hạn như các mặt tường với

vận tốc bằng không, các mặt đối xứng, vận tốc lối vào và áp suất lối ra của dòng

không khí.

Trong  nghiên  cứu  này,  vận  tốc dòng  khí  lối  vào  được  thay đổi  trong  mô

phỏng. Vận tốc được thay đổi trong khoảng 0 đến 50 m/s với từng nấc 5 m/s, nó

phù hợp với phạm vi thử nghiệm của các UAV bay trong điều kiện tốc độ thấp.

Những mô phỏng này được lặp đi lặp lại ở các góc tới là 0 đến 20 độ. Sau đó,

lực khí động được đo trong từng mô phỏng, để xác định hệ số của lực nâng và

lực cản, và được so sánh với kết quả lý thuyết.

Các kết quả mô phỏng đã được khảo sát trong một số gói phần mềm khác

nhau của  ANSYS.  ANSYS Fluent có thể  cung cấp một số hình đồ họa, chẳng

hạn như các phân bố áp suất và vận tốc. Mặt khác, ANSYS Structural cho phép

xác định các lực, chuyển vị, ứng suất và biến dạng của cánh.

Hình 20 cho thấy biểu đồ miền phân bố áp suất trong luồng không khí, khi

vận tốc lối vào được áp dụng 25m/s.

Hình 20. Biểu đồ phân bố áp suất  42

Như có thể thấy, các vùng áp suất cao xuất hiện ở mép trước và mặt dưới

của cánh. Bên cạnh đó, các khu vực áp suất thấp hơn xảy ra ở bề mặt trên của

cánh. Phân tích này là chính xác với các lý thuyết về hình thành lực nâng.

Tốc độ dòng khí cũng là một đặc tính quan trọng liên quan đến áp suất khí động.

Các miền phân bố tốc độ được thể hiện trong hình 21.

Hình 21. Biểu đồ phân bố vận tốc

Ở mép trước và sát bề mặt của cánh máy bay, vận tốc của dòng khí là gần

như bằng không. Tuy nhiên, các dòng khí tăng tốc thay đổi rõ ràng ở bề mặt trên

của cánh.

Hình  22  thể  hiện  ứng  suất  tương  đương  trên  cánh  sử  dụng  ANSYS

Structural.

Hình 22: Ứng suất tương đương trên cánh

Ứng suất đạt tối đa trên cánh ở phần gắn cố định vào thân máy bay. Ngoài

ra, các lực nâng và lực cản đã được xác định bằng tương tác kết cấu chất lỏng

43

(FSI) (Jin vcs 2002) [45]. Các thành phần lực tương ứng với vận tốc lối vào đã

được thu thập.  Từ  những dữ liệu này,  hai đồ thị của các mối quan  hệ giữa lực

nâng,  lực  cản  so  với  vận  tốc  tương  đối  giữa  các cánh  và  luồng  không  khí  lần

)

)

N

N

(   n ả c   c ự L

(   g n â n   c ự L

Kết quả mô phỏng Kết quả lý thuyết

Kết quả mô phỏng Kết quả lý thuyết

Vận tốc tương đối (m/s)

Vận tốc tương đối (m/s)

lượt được thể hiện trong hình 23, và hình 24.

Hình 23. Quan hệ giữa lực nâng  với vận tốc tương đối.

Hình 24. Quan hệ giữa lực cản  với vận tốc tương đối.

Các kết quả được so sánh với công thức lý thuyết bởi phương trình (46) và

(47), với góc tới vô cùng nhỏ được chọn (a≈0.025). Kết quả so sánh cho thấy sự

phù hợp tốt giữa lý thuyết và mô phỏng. Do đó các phương pháp phân tích đề

xuất đã chứng minh một sự thay thế khả thi để có được các lực khí động và các

hệ số bằng cách thao tác các kết quả từ mô phỏng ANSYS. Tuy nhiên, nghiên

cứu thêm cần được đề xuất để giảm sự khác biệt trong kết quả ở các điều kiện

nhất định, và để cho phép tính toán ma sát liên quan đến lực nâng và cản. Các hệ

0.156

0.056

LC 

DC 

quả cũng cho phép tính được rằng lực nâng gấp lực cản khoảng 22,5 lần.

Không chỉ dự đoán được các đặc tính khí động, mô hình cánh máy bay có

thể là một lựa chọn rất tốt cho sản xuất UAV thử nghiệm trong tương lai. Các dữ

liệu thu được có thể hỗ trợ trong các nghiên cứu trong tương lai chẳng hạn như

việc lựa chọn các cảm biến và thiết kế hệ thống điều khiển thích hợp.

số lực nâng và lực cản được xác định lần lượt là , . Từ kết

44

2.5.2. Tối ưu hình dạng khí động sử dụng phương pháp SQP

Trong những năm gần đây, sự phát triển của các máy bay không người lái

(UAV) đã dẫn tới nhiều nghiên cứu cho các ứng dụng cả quân sự và dân sự như

giám sát ven biển, quan sát thời tiết, theo dõi cháy rừng, thu thập dữ liệu khoa

học (Gupta vcs 2013[39],  Liu vcs 2014 [61],  Abdelkader vcs 2013 [9], Hardin

vcs 2011 [40], Rathinam vcs 2008 [73], Turner vcs 2012 [84]). Ở Việt Nam, lĩnh

vực này đã thu hút được sự chú ý đặc biệt bởi các nhà thiết kế và các nhà khoa

học, với sự hỗ trợ của các công cụ hiện đại. Tuy nhiên, hầu hết UAV bay ở một  số Reynolds thấp, trong khoảng từ 105 đến 106 (Thomas vcs 2003) [83]. Vì vậy,

tối ưu  hóa hình dạng thường rất khó  áp dụng bởi vì các lực  nhớt tác động lớn

hơn trên khí động của các thiết diện cánh tại số Reynolds. Việc tối ưu hóa hình

dạng khí động đã trở thành một phần thiết yếu của một công cụ tối ưu hóa thiết

kế cho các ứng dụng hàng không.

Từ ba thập kỷ qua, sự gia tăng hiệu quả tính toán và các tiến bộ mới trong

lĩnh vực CFD đã cho phép các nhà thiết kế mô hình và phân tích các hệ thống

phức tạp với một lượng thời gian tính toán hợp lý. Như kết quả của những tiến

bộ này, giờ đây ta có thể để kết hợp một lời giải CFD với một kỹ thuật tối ưu

hóa số  cho  thiết  kế  tối  ưu  của  các  hình  dạng  cánh  UAV.  Trong  lĩnh  vực này,

nghiên cứu đã tập trung vào việc giảm số lượng thời gian tính toán cần thiết để

tính các hàm chức năng và các gradient.

thường liên quan đến việc giải quyết một bài toán tối ưu hóa phi tuyến ràng buộc

(Venter 2010) [85]. Phương pháp  tuần tự tiểu cục có thể được sử dụng để giải

quyết các vấn đề phi tuyến ràng buộc cho một bài toán việc tối ưu phi tuyến ràng

buộc đơn giản. Chú ý rằng những phương pháp này là phương pháp tiểu cục, do

đó, các công cụ thiết kế không đảm bảo rằng hình dạng tối ưu là một tối ưu toàn

cục.  Các  phương  pháp  tuần  tự  tiểu  cục  thường được  sử  dụng là:  the  modified

method  of  feasible  directions  (MMFD)  (Thomas  vcs  2003)  [83],  sequential

linear  programming  (SLP),  sequential  quadratic  programming  (SQP),  and

Cũng giống như các ứng dụng kỹ thuật khác, tối ưu hóa hình dạng khí động

45

response  surface  approximation  methods  (RSM). Trong  đó,  phương pháp  SQP

được đánh giá là rất mạnh, có tốc độ hội tụ nhanh và hiệu quả cao nhất trong các

bài toán tối ưu phi tuyến ràng buộc này.

Để một  máy  bay  có  thể  bay được,  lực  nâng  khí  động  trên  các  cánh  phải

thắng được trọng lượng cơ thể máy bay, và lực cản phải nhỏ hơn lực đẩy của các

động cơ máy bay. Các lực khí động này thường phụ thuộc vào hình dạng thiết

diện  cánh,  vận  tốc  và  mật  độ  dòng  khí,  cấu  trúc  cánh  và  góc  nâng  khí  động.

Thông  thường,  trong  phần  lớn  các  máy  bay,  góc  nâng  tăng  dẫn  đến  lực  nâng

tăng. Tuy nhiên, góc này giới hạn bởi một góc tới hạn, mà tại đó lực nâng không

thể tăng được nữa mà giảm xuống đột ngột. Với các dòng khí tại số Reynolds

thấp, góc nâng tới hạn thường nhỏ.

Thiết diện Eppler 66, một loại cũng thường được sử dụng trong chế tạo các

cánh UAV ở số Reynolds thấp, sẽ được lựa chọn trong bài toán tối ưu này. Theo

các tài liệu mở về loại cánh này, ở số Reynolds thấp, khi góc nâng tăng từ 0 đến

5 độ, thì hệ số lực nâng tăng. Tuy nhiên, hệ số lực cản cũng tăng nhanh chóng,

dẫn tới hiệu quả thấp của hệ số lực nâng trên lực cản. Trong mục này, phương

pháp  SQP  cho  tối  ưu  hình  dạng  khí  động  trên  cánh  đơn  2  chiều  sẽ  được  thực

hiện. Bên cạnh đó, phần mềm ANSYS sẽ được sử dụng để tính toán và dự đoán

các đặc tính khí động trên các thiết diện cánh. Cuối cùng, nghiên cứu sẽ chỉ ra

và chứng minh một thiết diện cánh thỏa mãn các đặc tính khí động yêu cầu.

Tương tự  phần lớn các ứng dụng  kỹ thuật,  bài toán tối ưu  phi tuyến ràng

buộc được thực hiện nhằm mục đích:

Tối thiểu

trong khi

(48)

( ) 0

j



1,2,...,

q .

f x ( )



(

,

,...,

f x x 1 2

x )n

jg x   với

trong  đó,

f x   là  hàm mục  tiêu,  x là một vectơ  của các biến  thiết  kế, n là số

( )

lượng các biến thiết kế, và

( )

jg x  là các ràng buộc. Một chú ý rằng bài toán tối

đa  hóa  một  hàm

( )m x ,  sẽ  tương  tự  như  bài  toán  tối  thiểu  hóa  hàm

. Vì vậy, các bài toán này có thể giải bằng cách sử dụng thuật toán

f x ( )

 

m x ( )

Phương pháp tối ưu SQP

46

f x g x  là phi tuyến, liên tục và liên tục

( ),

( )

tối ưu tương tự. Giả sử rằng các hàm

R



)



0} ,

*x  là một kết quả tối

với  các  vi  phân  bậc  nhất,  bậc  hai.  Sau  đó,  không  gian  biến  được  xác  định  là

g x   ( {x : j

x R  được coi là một điểm khả thi. Và

*( f x

)



f x

( ),

    x R .

thiểu của bài toán tối ưu nếu

Thuật  toán lập trình toàn  phương  liên tiếp  Sequential  Quadratic

Programming (SQP) được coi là một trong những phương pháp hiệu quả nhất để

giải quyết bài toán tối ưu hóa ràng buộc phi tuyến.  Để thực hiện phương pháp

này, yêu cầu chính là một lời giải lập trình bậc hai hiệu quả. Giải pháp tối ưu của

bài toán phải đáp ứng các điều kiện Karush-Kuhn-Tucker (KKT) [Antoniou vcs

*

j ) 0  for





1,2,...,

q ;

g x ( j

q

*

*

*

(

x

* 

,

)

 

f x (

)



 ) 0;

  x

x

*   j

x

g x ( j



j

 1

2003] [18]:

*

j ) 0 for





1,2,...,

q ;

*   j

g x ( j

x



j 0  for



1,2,...,

q ;

*  j

q

*

*

*

(49)



(

x

* 

,

)



f x (

)

)

(

* * x  ,

)

*  j

g x ( j

  x

 

j

 1

*

Trong  đó là  một  hàm  Lagrange, là

(

* x   là cặp biến tối ưu.

)

,

gradient của hàm Lagrange và

Nó là cần thiết để sử dụng các phương trình đơn giản hơn và sau đó giải quyết

điều kiện KKT, mà có thể được viết dưới dạng ma trận như sau:

Tối thiểu

trong khi

(50)

 ;  0

H

p

b  A k x k

T    k x

T x

x

k

1 2

q

với

H

 

)



)

k

2 x

f x ( k

k

j

2 x

g ( j

x k

x   là  số  gia  của  các  biến  thiết  kế

;kx

 u ( )

j

 1



)

)

T x

được gọi là ma trận Hessian;

và

 

);

;

.

p k

f x ( k

x

A k

b k

)

)

T x

g x 1( k ... g x ( q k

g x 1( k ... g x ( q k

     

    

     

    

các vấn đề một cách tuần tự. Một tập mới của các phương trình tương đương với

47

0;1a 

x  với

Tham số α được sử dụng như kích thước vec tơ . Bài toán được

giải với α ban đầu, vec tơ các biến được cập nhật sử dụng biểu thức:

x 1k

 

 x a k k

(51)

 1

Và nhân tử Lagrange cũng được cập nhật bằng cách sử dụng biểu thức:





 k

 

1

T A A k k

A k

 a  H k x

p k





(52)

Tại  lần  lặp  đầu  tiên,  ma  trận  Hessian  của  hàm  Lagrange  được  xấp  xỉ  bởi  các

dạng ma trận ban đầu. Trong lần lặp tiếp theo, tính xấp xỉ của ma trận Hessian

kH  1,

thường được sử  dụng trong phương pháp quasi-Newton,  được tính bằng

k

công thức  BFGS:

H

H





 

1

k

k

T H d d H k k k d H d

T k

k

k

T  k k T  d k k

(53)

   y

;

d







;

   1

 

k

a x

x k

 1

x k

k

k

H d k

k

y

,

,

;









k

   x

x k

 k

 1

 1

   x

x k

 k

 1



0.2

d H d

T d y k k

k

k

T k

trong đó,



if



0.2

d H d

T d y k k

k

k

T k

d H d k k T  d y k k

T 0.8 k T d H d k

k

k

 1                        if       

.

Tại đây, phương pháp SQP thực hiện trong gói các công cụ thiết kế tối ưu

hóa (DOT), là một trong những phiên bản phổ biến nhất của SQP. Phương pháp

I  với  I  là ma trận ban đầu,

kx  và

kH ,

Bước 1: đặt  Bước 2: tính  f x

x 0  g x

0, k    k

và  k

Bước 3: tính gradient của các hàm mục tiêu và các ràng buộc

Bước 4: giải phương trình (49),

Bước  5:  tính  các  biến  thiết  kế  sử  dụng  phương  trình  (51)  và  các  nhân  tử

Lagrange sử dụng phương trình (52),

Bước 6: cập nhật xấp xỉ ma trận Hessian sử dụng phương trình (53),

này có thể được tóm tắt trong các bước sau đây:

48



 f x k



0,001

 f x k 

   1  f x k

Bước 7: kiểm tra độ hội tụ. Nếu trong hai lần lặp liên



  i

0,001;

1, 2,...,

n

 f x 

1k

id sẽ cần được tính toán,  Bước 8: cập nhật vòng lặp

k

k   và quay trở lại bước 2.

1,

tiếp hoặc thì hội tụ đặt được, và STOP. Chú ý

Trong các phần tiếp theo, phương pháp SQP sẽ được áp dụng cho bài toán

tối ưu hóa hình dạng cánh máy bay.

Một máy bay để bay được cần thỏa mãn các yêu cầu:

Nâng ≥ Trọng lượng

Cản ≤ Đẩy

Vì vậy, một cánh mong muốn là có thể đáp ứng một yêu cầu lực nâng nhất

định, và có một lực cản tối thiểu tại cùng một thời điểm. Từ phương trình (50),

một bài toán tối ưu hóa hình dạng khí động được đưa ra:

  và  0

( )

  c x l

desired c l

jg x  ;  ( ) 0

dc x  trong khi

Tối thiểu



 x  lần lượt là các hệ

  lc x  và

dc

desired

với,  x  đại diện cho hình dạng thiêt diện cánh,

lc

số lực nâng và cản của thiết diện cánh, là hệ số lực nâng tối thiểu yêu cầu

 

jg x  là các ràng buộc hình học của thiết diện cánh. Trong nghiên cứu này,

và

các  hàm  mục  tiêu  và  các  ràng  buộc  khí  động  thu  được  bằng  cách  giải  các

phương trình RANS (compressible mass weighted Reynolds Averaged Navier-

4

5.10 ,

Đặc biệt, các ràng buộc hình học được áp

Stokes) của dòng khí xung quanh cánh. Bài toán thiết kế là để có được một cánh

5 . a 

đặt  để  có được  cánh  độ  dày  tối  thiểu  là  1%  của  chiều  dài dây  cung  thiết diện

cánh.

với lực cản tối thiểu và hệ số nâng tối thiểu là 0,8, tại một số Reynolds là   và với một góc tới là 5 độ 

49

Dạng ban đầu của thiết diện (Eppler 66)

Dạng thiết diện tối ưu (SQP)

Hình 25. Hình dạng thiết diện cánh máy bay ban đầu và cánh tối ưu

Để giải bài toán thiết kế, thiết diện cánh Eppler 66 được sử dụng như thiết

diện ban đầu. Thiết diện này được lựa chọn bởi vì nó là một trong các thiết diện

được nhắc tới trong một số thiết kế cánh máy bay tại số Reynolds thấp. Kết quả

lời giải tối ưu có được sử dụng thuật toán SQP được thể hiện trên hình 25. Trên

đó, hình dạng cánh ban đầu và cánh tối ưu được đưa ra. Qua đó có thể thấy rằng

để có được cánh tối ưu với lực cản tối thiểu, thiết diện cánh cần được làm mỏng

hơn.

Trong nghiên cứu này, ANSYS workbench 16.0 được sử dụng để thiết kế

thước nàu được sử dụng để chế tạo mô hình cánh thực nghiệm, nó cũng phù hợp

với các dữ liệu mở của một số mẫu UAV đã chế tạo tại Việt Nam. Để chạy mô

phỏng, các điều kiện biên và phương trình giải được cài đặt tương tự như các mô

phỏng trên, với vận tốc lối vào là 25 m/s và góc tới tại

5 ,a   lần lượt cho hai

mô hình cánh ban đầu và tối ưu.

Các hệ số lực nâng và lực cản lần lượt cho các cánh máy bay ban đầu và

cánh máy bay tối ưu được tính toán được thể hiện trong bảng 2. Như có thể thấy

một tập hợp các đặc điểm của một thiết diện cánh UAV. Các thông số được lựa  chọn,  như  chiều  dài  dây  cung  là  0,4  m  và  diện  tích  cánh  là  1,2  m2.  Các  kích

50

rằng sự ràng buộc lực nâng là hiệu quả, bởi vì các hệ số lực nâng vẫn còn lớn

hơn giá trị 0,8. Phương pháp SQP giảm đã giảm lực cản tới gần 20% so với cánh

máy bay ban đầu với. Do đó, tỷ lệ lực nâng trên lực cản tăng từ 15,61 đến 19,05.

Bảng 2: Các đặc tính khí động lực trên các thiết diện cánh tại Re=5.104 và a=5o

Đặc tính Eppler 66 SQP

Hệ số lực nâng 0,826 0,804

Hệ số lực cản 0,0529 0,0422

Hệ số lực nâng trên lực cản 15,61 19,05

Phân tích CFD được thực hiện sử dụng ANSYS Fluent. Mô hình rối  k 

được  sử  dụng  để  phân  tích.  Sự  phân  bố  áp  suất  và  vận  tốc  của  các  tác  động

không khí trên cánh máy bay, được thể hiện trong các hình dưới đây.

Hình 26.  Hình ảnh phân bố áp suất tại Re=5.104 và a=5o trên thiết diện cánh  Eppler 66  (hình trái) và thiết diện tối ưu theo SQP (hình phải).

Hình  26  cho  thấy sự phân bố áp suất khí động trong các  trường dòng khí

cách sử dụng phương pháp SQP. Đối với thiết kế tối ưu, có thể thấy rằng áp suất

tối đa tại mép trước của cánh giảm, từ 364 Pa xuống còn 282 Pa. Việc giảm áp

suất tối đa tại mép trước là lý do chính cho thấy  việc giảm lực cản. Ngoài ra, áp

suất lên bề mặt dưới của cánh là cao hơn, trong khi áp lực trên bề mặt phía trên

là nhỏ hơn so với cánh máy bay ban đầu. Độ chênh lệch áp suất giữa mặt dưới

và mặt trên vẫn đảm bảo yêu cầu hình thành lực nâng. Tương tự, phân bố vận

tốc được  thể  hiện  trên hình  27.  Thang  màu  sắc  từ  màu xanh  đậm đến màu đỏ

đậm thể hiện  độ lớn của vận tốc.  Trên mép trước và sát bề mặt của cánh máy  51

xung quanh cả hai cánh máy bay ban đầu và cánh máy bay tối ưu thu được bằng

bay, vận tốc của dòng khí là gần như bằng không. Ngược lại với sự phân bố vận

tốc ở bề mặt trên của cánh máy bay tối ưu là cao hơn, trong khi đó vận tốc ở bề

mặt dưới là nhỏ hơn so với cánh máy bay ban đầu.

Hình 27.  Hình ảnh phân bố vận tốc Re=5.104 và a=5o trên thiết diện cánh  Eppler 66 (hình  trái) và thiết diện tối ưu theo SQP (hình phải).

Hình 28 cho thấy ứng suất tương đương trên các cánh máy bay, có được từ

sử dụng gói ANSYS Structural. Hợp kim nhôm 7075 T6 cũng đã được sử dụng

trong  các  thiết  kế  của  nghiên  cứu  này.  Các  đặc  tính  vật  liệu  liên  quan  có  thể

được tìm thấy trong rất nhiều cuốn sách cơ bản. Như có thể thấy rằng ứng suất

đạt tối đa ở phần gắn vào thân máy bay, giảm xuống từ cánh ban đầu đến cánh

tối ưu. Hơn nữa, sự phân bố ứng suất của cánh tối ưu là đồng đều so với cái ban

đầu. Như đã biết ứng suất liên quan đến độ bền của cánh.

Hình 28: Hình ảnh phân bố ứng suất trên cánh tại Re=5.104 và a=5o đối với thiết diện cánh   Eppler 66 (hình trái) và thiết diện tối ưu theo SQP (hình phải).

Các thuật toán SQP đã thể hiện được là một thuật toán tối ưu hóa hiệu quả

nhất do  thực tế rằng nó  yêu cầu số lượng thấp nhất về các tính toán hàm  chức

năng. Trong trường hợp bài toán tối ưu hình dạng cánh máy bay này, nó chỉ đòi  52

hỏi có 26 tính toán chức năng. Cánh máy bay tối ưu có được đã chứng minh có

các đặc tính hiệu suất tốt hơn bằng cách sử dụng gói công cụ ANSYS.

Như  vậy  trong  mục  này,  một phương  pháp để  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.

2.5.3. Mô phỏng CFD trên cánh máy bay với các góc tới lớn

Mục  này  mô  tả  quá  trình  hoàn  thành  mô  hình  hóa  và  mô phỏng  của  các

luồng không khí xung quanh cánh máy bay. Dựa trên lý thuyết khí động lực, các

phân tích CFD thu được bằng cách sử dụng gói phần mềm ANSYS Fluent. Các

bộ giải dựa trên thuật toán SIMPLE đã được lựa chọn. Các mô phỏng đã được

tiến hành trên mô hình cánh có thiết diên cánh là NACA 2412. Các tính toán khí

động  lực  được  thực  hiện  đối  với  các  góc  tới  khác  nhau.  Mục  đích  chính  của

nghiên cứu trong mục này này để dự đoán được các đặc tính khí động điển hình

được góc tới giới hạn, mà tại đó lực nâng của cánh đạt hiệu quả lớn nhất.

Động lực học chất lỏng tính toán (CFD) [Petrila vcs 2005] [66] là phân tích

hệ  thống  liên  quan  đến  dòng  khí  chất  lỏng  khí  bằng  cách  thực  hiện  các  mô

phỏng dựa trên máy tính. Ngày nay, CFD được ứng dụng ngày càng nhiều trong

lĩnh vực hàng không, từ thiết kế, chế tạo cho đến nghiên cứu và phát triển máy

bay.  Máy bay  như  một  đối tượng  có thể bay trong  các  luồng  chảy  không  khí.

Điều này có được là nhờ các lực nâng khí động đã thắng trọng lượng cơ thể của

máy bay.

và cải thiện hiệu suất cánh máy bay ở góc tới lớn. Kết quả tính toán đã xác định

53

Các lực khí động hình thành do chuyển động tương đối của vật thể bay với

dòng  khí.  Hai  thành  phần cơ  bản gồm  lực  nâng  vuông  góc  với  hướng  chuyển

động, và lực cản song song với hưởng chuyển động [Mohamed 2000] [64]. Các

lực này được đặc trưng bởi hệ số lực nâng và hệ số lực cản tương ứng. Các hệ số

khí động này chủ yếu phụ thuộc hình dạng cánh, thiết diện cánh tham chiếu, mật

độ, độ nhớt và vận tốc của dòng không khí và góc tới của cánh máy bay.

Góc tới thường được định nghĩa là góc giữa dây cung thiết diện và hướng

chuyển động tương đối của cánh máy bay với dòng khí. Phần lớn các cánh yêu

cầu  góc tới dương  để phát  sinh lực  nâng  [Rumsey  vcs  2002]  [75]. Nói chung,

khi tăng góc tới thì hệ số khí động tăng. Tuy nhiên, góc tới này là có giới hạn,

nếu vượt qua góc tới giới hạn hệ số lực nâng giảm, trong khi hệ số lực cản tăng

nhanh chóng, gây nguy hiểm cho máy bay. Nghiên cứu này được thực hiện để

xem xét hiệu quả khí động khi tăng góc tới và dự  đoán góc tới giới hạn trong

trường hợp cụ thể, sử dụng công cụ mô phỏng và tính toán CFD.

Trước hết ta trình bày lại một số các kết quả lý thuyết.

Một cơ hệ đơn giản của mô hình cánh máy bay cứng, thiết diện cánh không

đổi.  Mặt  phẳng  thiết  diện  được  treo  bằng  một  lò  xo  xoắn  gắn  với  tường  của

Biên dẫn

đường hầm gió, được thể hiện như hình 29 [Dowell vcs 2015] [31].

O

L

e

α

U

k

MAC

Biên theo

W

ab

c

Hình 29. Phác họa thiết diện đặc trưng của cánh máy bay.

54

Độ  cứng  xoắn  của  lò  xo  được  giả  định  là  k.  Giả  sử  vật  thể  có  thể  xoay

quanh một điểm  nối với lò  xo,  được  gọi  là  tâm đàn hồi  của  cánh. Góc  tới  khí

động tổng cộng của cánh là a, nó bằng tổng của góc tới ban đầu a0 (khi lò xo

không xoắn), cộng thêm một lượng do sự xoắn của lò xo ae.





a a a 0 e

(54)

Một điểm khác trên cánh máy bay gọi là tâm khí động, tại đó mômen khí

động là độc lập với góc tới. Vì vậy, phương trình cân bằng mômen của cơ hệ đối

với tâm đàn hồi được cho là:



 Le Wd

 M M y

AC

(55)

trong đó,

My=kae: mômen đối với tâm đàn hồi

MAC: mômen đối với tâm khí động, cả hai mômen này được quy ước là chiều

dương quay cùng chiều kim đồng hồ.

L: lực nâng, quy ước chiều dương hướng lên

W: trọng lượng cánh

e,d: lần lượt là khoảng cánh từ tâm khí động, trọng tâm cánh đến tâm đàn hồi.

L C qS



Từ lý thuyết khí động [Anderson 2010] [12], ta có:

L 

M

AC

C qSc MAC

(56)

C



C



a

L

L 0

C  L  a

: là một hệ số mômen tâm khí động, được coi không đổi.

C

C

MAC

MAC 0

Chú ý rằng, các hàm vô hướng

được chọn tại góc tới

0a .

,LC

0

0MACC

2

gọi là áp suất khí động, với mật độ dòng khí , và vận tốc dòng khí

q



 U 2

tương đối U.

c chiều dài dây cung thiết diện, l sải cánh máy bay, và diện tích cánh tham chiếu

được xấp xỉ  S

  .  c l

trong đó là hệ số lực nâng

55



MAC

Giải phương trình (56) tìm độ uốn đàn hồi, ta được:



a e

 

a C qSc C qS e Wd L 0 k C qSe L

(57)

Từ phương trình (57), ta thấy rằng tâm khí động ở phía trước tâm đàn hồi

khi e>0, làm cho kết quả tiến tới giới hạn khi mẫu thức tiến tới vô cùng nhỏ hay

ae tăng lên khi q tăng đáng kể. Mẫu của biểu thức tính ae là một loại độ cứng tác

dụng, nó giảm khi q tăng. Khi mẫu thức triệt tiêu, hiện tượng phân kỳ xảy ra, khi

đó áp suất khí động được xác định:

q



D

k C Se L

(58)

Như vậy, rõ ràng là khi tâm khí động trùng với tâm đàn hồi, nghĩa là e=0,

áp suất khí động phân kỳ tiến tới vô hạn. Còn nếu tâm khí động lực ở phía đuôi

sau tâm đàn hồi, nghĩa là e<0, thì áp suất khí động phân kỳ trở thành âm. Do các

yêu cầu vật lý, áp suất khí động phải dương và hữu hạn, nên trong trường hợp

phân kỳ này là không thể.

C

C

  và  0

d  , ta được:

0

L 0

MAC 0

e

a 0

Để đơn giản, ta giả sử rằng



a e

qS k

q

1



 C L  a Se C  L  a k

(59)

k



q Se D

C  L  a

có thể được đơn giản như sau:

(60)

a  e

 1

q q /  

 a D 0  q q /

D

Ta có lực nâng tỉ lệ với góc tới. Vì vậy, độ thay đổi lực nâng được cho bởi:

(61)





 

1

 / q q D  q q /



 L L 0

a e a 0

D

Kết hợp phương trình (59) với (58), ta có thể cho ; do đó, ae

56



q

q . Đối với

L 0

D

LC  qS a 0  a

Với . Cả ae và L đều tiến tới vô cùng khi



 2

LC   a

tấm cứng trong dòng khí hai chiều không nén được, ta có .

Tóm lại, từ phương trình (61) cho thấy lực nâng tăng khi góc tới tăng. Tuy

nhiên,  lực  nâng  chỉ  tăng  tới  giới  hạn  khi  áp  suất  khí  động  tăng  tới  giới  hạn.

Ngoài ra, hệ số lực nâng không phải là không thay đổi, nó phụ thuộc nhiều yếu

tố  môi  trường  xung  quanh.  Bên  cạnh  đó,  cấu  trúc  cánh  cũng  không  thể  chịu

được biến dạng vô hạn, nó dễ xảy ra hư hỏng cấu trúc khi góc tới giới hạn. Điều

đó rất nguy hiểm nên áp suất khí động thường giới hạn dưới mức phân kỳ. Các

phần tiếp theo sẽ thực hiện các mô phỏng trên cánh  máy bay, để tính toán  lực

nâng và lực cản khi thay đổi góc tới, sử dụng các gói phần mềm ANSYS.  Công

cụ mô phỏng này cho phép dự đoán được các đặc tính khí động có độ chính xác

tốt  và  hiệu  quả  nhanh.  Bây  giờ  ta  trình  bày  các  kết  quả  mô  hình  hóa  và  mô

phỏng. Thiết diện cánh NACA 2412 [Eppler 1990] được lựa chọn để tạo ra mô

hình  cánh  máy  bay,  sử  dụng  các  phần  mềm  CAD  như  AutoCAD,  SolidWork,

CATIA  geometry,  có  thể nhập vào  ANSYS Geometry. Mô  hình  cấu  trúc  cánh

tạo thành được thể hiện trên hình 30. Mô hình sau đó được tạo lưới bằng các gói

tích hợp sẵn trong ANSYS.

Hình 30. Hình học của hình dạng cánh

Giả sử cánh được làm bằng hợp kim nhôm 7075 T6, bay trong dòng không

khí. Các thông số liên quan đến vật liệu và tính chất cánh đã có trong nhiều các

tài liệu mở. Để chạy mô phỏng, các điều kiện biên và thông số ban đầu được cài

đặt.  Vận  tốc  dòng  khí  tương  đối  lối  vào  được  chọn  là  240  m/s,  tức  là  ở  số

57

Reynolds  lớn;  trong  khi  góc  tới  được  thay  đổi  0  đến  20  độ,  ứng  với  các  mô

phỏng khác nhau.

Phân bố áp suất (Pa) Phân bố vận tốc (m/s)

Hình 31. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh tại góc tới  là 0 độ.

Phân bố vận tốc (m/s) Phân bố áp suất (Pa)

Hình 32. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh tại góc tới  là 10 độ.

Phân bố áp suất (Pa)

Phân bố vận tốc (m/s)

Hình 33. Biểu đồ phân bố áp suất (trái) và phân bố vận tốc dòng khí xung quanh cánh tại góc tới  là 18 độ.

58

Phân tích CFD được thực hiện sử dụng ANSYS Fluent. Mô hình chảy rối

k- được sử dụng trong phân tích. Phân bố áp suất khí động và vận tốc dòng khí

tác động xung quanh cánh được rút ra trong một số trường hợp khi thay đổi góc

tới, và được thể hiện trong các hình bên trên từ hình 31 đến hình 33.

Độ lớn áp suất hay vận tốc được thể hiện bằng các miền màu sắc từ xanh

đậm đến đỏ đậm, tương ứng với giá trị từ nhỏ đến lớn. Các kết quả mô phỏng

cho thấy, áp suất ở mặt trên là thấp hơn áp suất mặt dưới, tạo ra độ chênh lệch

áp suất, từ đó làm lực nâng hướng lên trên, phù hợp với lý thuyết hình thành lực

nâng. Ngược lại, vận tốc dòng khí phía mặt trên lại cao hơn so với mặt dưới. Khi

tăng góc tới, nói chung độ lớn áp suất và vận tốc khí động đều tăng, các miền

màu thể hiện càng rõ rệt. Độ chênh lệch áp suất tăng, dẫn tới lực nâng tăng. Tuy

nhiên, ở các góc tới lớn, dòng khí có hiện tượng chuyển sang trạng thái hỗn độn,

điều đó  làm lực  nâng  và  có thể dẫn  tới  nguy  hiểm. Hơn  nữa,  áp  suất  lớn nhất

thường  là  ở  mép trước  của  cánh,  nó  thể  hiện  thành  phần  lực cản tác  động  lên

cánh.  Khi  tăng  góc  tới  lớn,  áp  suất  này  tăng  theo  rất  nhanh,  cho  thấy  lực  cản

cũng tăng mạnh, càng làm giảm hiệu quả bay.

Các mô phỏng để tính toán lực nâng và lực cản khi góc tới là 0 độ, 5 độ, 10

độ và từ 12 đến 20 độ cũng được thực hiện. Kết quả được đưa vào bảng 3.

Bảng 3. Các giá trị lực nâng và lực cản tương ứng với các góc tới khác nhau  TT Góc tới AoA (độ) Lực nâng L (N) Lực cản D (N) 2.7378e5  1  1.1446e6  2  2.0844e6  3  2.3521e6  4  2.4916e6  5  2.6377e6  6  2.7481e6  7  2.8487e6  8  2.9152e6  9  2.9795e6  10  2.8981e6  11  2.8022e6  12

0.9223e4  1.1254e4  1.3231e4  1.8965e4  2.2818e4  2.6742e4  3.2974e4  4.5845e4  7.4378e4  1.8436e5  3.6136e5  6.4626e5 0  5  10  12  13  14  15  16  17  18  19  20

59

Ta thấy rằng, khi góc tới tăng, lực nâng tăng song không phải tăng mãi mãi.

Bởi sau khi qua góc tới 18 độ, lực nâng bắt đầu giảm,  ta có thể nhận biết đây

chính là góc tới giới hạn. Tuy nhiên, với lực cản, lực này tăng nhẹ ở góc tới nhỏ,

nhưng sau đó tăng lên nhanh chóng khi góc tới lớn, và tiếp tục tăng  mạnh khi

góc tới vượt qua giới hạn.

Các hệ số khí động được tính toán và thể hiện lần lượt trên hình 34 và hình

35. Các kết quả này được so sánh với các tính toán trước đó của NACA, với số

Reynolds khoảng 10e6. Có thể thấy rằng, kết quả mô phỏng là phù hợp. Góc tới

giới hạn được xác định vào khoảng 18 độ, tại đó hiệu quả lực nâng có thể là tốt

Kết quả mô phỏng

Kết quả lý thuyết

g n â n   c ự l   ố s   ỉ

T

Góc tới (độ)

nhất.

Hình 34. Đồ thị hệ số lực nâng theo góc tới

Kết quả mô phỏng

Kết quả lý thuyết

n ả c   c ự l   ố s   ỉ

T

Góc tới (độ)

Hình 35. Đồ thị hệ số lực cản theo góc tới

60

Tóm  lại,  trong  mục này  đã trình  bày các kết  quả  nghiên  cứu này đã  thực

hiện phân  tích CFD  trên mô hình cánh  máy bay có thiết diện NACA 2412,  sử

dụng  các  gói  phần  mềm  ANSYS.  Các  mô  phỏng  đã  được  thực  hiện  tại  số

Reynolds lớn, và một số góc tới khác nhau. Các đặc tính khí động như phân bố

áp suất, vận tốc dòng khí và các lực khí động đã được tính toán và đưa ra. Từ

phân bố áp suất và vận tốc dòng khí cho thấy kết quả mô phỏng phù hợp với lý

thuyết khí động và hình thành lực nâng. Từ phân tích lực khí động cho thấy, lực

nâng hay hệ số lực nâng tăng lên khi tăng góc tới đến 18 độ, nhưng sau góc đó

thì chúng lại giảm. Đối với lực cản và hệ số lực cản, chúng tăng nhẹ ở góc tới

nhỏ và tăng mạnh ở các góc tới lớn. Các kết quả mô phỏng được so sánh với lý

thuyết  tính  toán  trước  đó, cho  thấy  sự tương  đồng  tốt.  Như vậy,  phương  pháp

mô  phỏng  số  này  có  thể  dự  đoán  hiệu  quả  các đặc  tính  khí  động  và  xác  định

được góc  tới  giới  hạn,  nếu  vượt  qua  góc tới giới hạn  này  có  thể  dẫn tới nguy

hiểm  cho  máy  bay.  Phương  pháp  mô phỏng  số  này  cũng  hứa  hẹn  sẽ  được  áp

dụng hiệu quả trong các nghiên cứu tiếp theo.

Kết luận chương 2

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 chuyển động trong dòng khí không nén

được. Phương trình phi tuyến thu được từ mô hình được dùng để phân tích đáp

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 và tính

toán số cũng như  phương pháp CFD và CSM để phân tích các hiện tượng khí

động, khí đàn hồi và dao động flutter xuất hiện trên thiết diện cánh.

ứng cũng như các hiện tượng flutter. Sau khi thiết lập phương trình dao động hai

61

CHƯƠNG 3. PHÁT TRIỂN KỸ THUẬT ĐỐI NGẪU CHO BÀI TOÁN

DAO ĐỘNG PHI TUYẾN

3.1. Phương pháp tuyến tính hóa tương đương

Các  hệ kỹ  thuật  thường là các  hệ  phi  tuyến  và  mô  hình tuyến tính  là  mô

hình  đơn  giản  hay  được  dùng  nhất  để  mô  tả  gần  đúng  các  hệ  này.Tuy  nhiên

trong nhiều trường hợp cần các mô tả chính xác hơn người ta vẫn phải sử dụng

mô hình phi tuyến và các phương pháp giải tích xấp xỉ là các công cụ thích hợp

để phân  tích  các bài toán phi tuyến. Ngày  nay, cùng  với  sự  phát  triển  của  các

phương pháp số, sự kết hợp với phương pháp giải tích giúp cho việc giải các bài

toán phi tuyến trở nên thực tế và khả thi hơn. Một trong những phương pháp giải

tích được sử dụng phổ biến là phương pháp tuyến tính hóa tương đương, do N.

Krylov  and  N.  Bogoliubov  (1937)  [48]  đề  xuất.  Ý  tưởng  cơ  bản  của  phương

pháp tuyến tính hóa là thay thế hệ phi tuyến ban đầu bằng một hệ tuyến tính mà

một  số  các  tính  chất  động  lực  học  của  hệ  phi  tuyến  có  thể  được  nghiên  cứu

thông  qua  hệ  tuyến  tính  này.  Do  vậy,  dẫn  tới  khái  niệm  về  hệ  tuyến  tính  hóa

tương đương.

Phương pháp tuyến tính hóa tương đương của Krylov và Bogoliubov đã được

nhiều  nhà  khoa  học  phát  triển  cho  lĩnh  vực  dao  động  phi  tuyến,  sau  đó  được

Caughey  phát  triển  trong  đó  tiêu  chuẩn  tương  đương  dựa  trên  điều  kiện  bình

1990 [74], Socha 2008 [77], Proppe vcs 2003 [70], Langley 1988 [49]). Ý tưởng

cơ bản là thay thế dao động phi tuyến ban đầu bằng dao động tuyến tính có cùng

kích động, trong đó các hệ số tuyến tính hóa được xác định bằng tiêu chuẩn cực

tiểu  hóa  trung  bình  bình  phương  của  sai  số  phương  trình.  Tiêu  chuẩn  do

Caughey  đề  xuất  thường  được  coi  là  cơ  sở  để  so  sánh  cho  các  phát  triển  lý

thuyết  của  phương  pháp tuyến tính  hóa tương  đương  về  sau. Vì  lí do  đó,  tiêu

chuẩn này còn được biết đến với các tên gọi khác như ‘tiêu chuẩn kinh điển’,

‘tiêu chuẩn thông thường’.

phương tối thiểu được sử dụng phổ biến nhất (Caughey 1963 [21], Roberts vcs

62

Như vậy, ý tưởng cơ bản của tiêu chuẩn tuyến tính hóa tương đương kinh

điển  cho  rằng  việc  cực  tiểu  hóa  sai  số  phương  trình  hay  sai  số  giữa  hàm  phi

tuyến và hàm tuyến tính tương đương có thể làm giảm thiểu sai số của nghiệm

xấp  xỉ  so  với  nghiệm  chính  xác.  Tuy  nhiên  cho  đến  nay,  vẫn  chưa  có  chứng

minh lý  thuyết về sự liên hệ sai số phương trình và sai số của nghiệm xấp xỉ. Có

nhiều  phát  triển  về  các  tiêu  chuẩn  tương  đương  khác  như  tiêu  chuẩn  tương

đương  năng  lượng  (Roberts  vcs  1990)  [74],  tiêu  chuẩn  tuyến  tính  điều  chỉnh

(Elishakoff  2009)  [33].  Các  kết  quả  cơ  bản  của  phương  pháp  tuyến  tính  hóa

tương đương được trình bày trong các sách chuyên khảo của Roberts và Spanos

(1990)  [74],  của  Socha  (2008)  [77],  đặc  biệt  bài  báo  tổng  quan  của  Crandall

(2006) [27] đã tổng kết 50 năm phát triển của phương pháp tuyến tính hóa tương

đương.

3.1.1. Tiêu chuẩn tương đương kinh điển

Để trình bày ý  tưởng cơ bản và một số phát triển của phương pháp tuyến

tính hóa tương đương, ta xét dao động phi tuyến được mô tả bởi phương trình



được biểu diễn dưới dạng

 x



2

 hx





,



f

 x g x x



  t

2  0

(62)

trong đó  x ,  x  và  x lần lượt là dịch chuyển, vận tốc và gia tốc của dao động;  h

( , )

g x x  là hàm phi tuyến đối với  x  và  x , f(t) là kích động

là hệ số giảm chấn,

cưỡng bức có thể là hàm tiền định hay ngẫu nhiên theo thời gian. Phương trình

(63)

 x



2

 hx



 x bx 



kx



f t ( )

2  0

trong đó b, k là các hệ số tuyến tính hóa. Sai số của phương trình sẽ là

(64)

                                              ( , ) e x x



 ( , ) g x x



 bx



kx

Theo tiêu chuẩn tương đương kinh điển ta yêu cầu sai số phương trình giữa

(62) và (63) phải thỏa mãn điều kiện cực tiểu của trung bình bình phương sai số

(65)

S



 ( , ) g x x



 bx



kx





2

k d

min b k ,

tuyến tính hóa tương đương của (62) có dạng

63

trong đó <(.)> là toán tử  trung bình  trên  một chu kỳ  khi nghiên cứu  dao động

tiền định, hoặc là toán tử trung bình theo xác suất khi nghiên cứu dao động ngẫu

nhiên. Từ (65) ta có



0;



0

 S kd  b

 S kd  k

(66)

xx 

0

,



Với giả thiết , giải hệ phương trình (66) thu được

b



   xg x x 2  x

,



(67)

k



  xg x x 2 x

(68)

Các hệ số tuyến tính hóa trong các công thức (67), (68) phụ thuộc vào các

đáp ứng chưa biết  x ,  x  là nghiệm của hệ tuyến tính tương đương (63). Ta thấy

 xuất hiện trong các công thức (35-36) cho thấy các

g

(x, x)

rằng hàm  phi tuyến

tính chất phi tuyến được thể hiện  trong việc xác định các  hệ số tuyến tính hóa

tương đương. Do vậy, có thể chờ đợi việc hệ tuyến tính hóa tương đương phản

ánh một số tính chất phi tuyến cơ bản của hệ phi tuyến gốc.

3.1.2. Tiêu chuẩn sai số thế năng

Elishakoff  và  Zhang  X  đã  đề  xuất  phương  pháp  tuyến  tính  hóa  tương



g x ( )

đương dựa trên cực tiểu  hóa  hàm  sai số thế năng cho trường hợp phương trình

1990) [74]. Hệ  số  tuyến  tính  hóa  k  trong phương  trình  (63) được xác định bởi

điều kiện cực tiểu sai số giữa thế năng U(x) của phần tử đàn hồi phi tuyến và thế

năng của phần tử  đàn hồi tuyến tính hóa tương đương

dưới dạng trung

kx

2 / 2

bình bình phương

2

2

(69)

U x ( )



kx



min k

1 2

  

  

trong đó

(62)  là  phương  trình  phi  tuyến  theo  dịch  chuyển,  g( , ) x x (Roberts  vcs

64

x

U x ( )

g

da a ( )

 

0

(70)

2

2

  2 x U x

Tiêu chuẩn (69) dẫn đến kết quả

k





  2 x U x 2

4

2

x

3

x

(71)

Tiêu chuẩn (69) áp dụng cho các hệ đàn hồi phi tuyến thu được kết quả khá tốt,

tuy nhiên không áp dụng được cho các hệ cản phi tuyến.

3.1.3 Tiêu chuẩn tương đương điều chỉnh

Anh  và  Di  Paola  đề  xuất  tiêu  chuẩn  tuyến  tính  hóa  tương  đương  chuẩn

Gauss có điều chỉnh dựa trên đánh giá: Sự thay thế hệ phi tuyến gốc bằng một

hệ tuyến tính  tự  nó đã làm  giảm mức độ phi tuyến (Elishakoff vcs 2009) [33].

Do  đó,  để đạt được nghiệm chính xác hơn thì cần cải thiện sự  giảm thiểu tính

phi tuyến nêu trên. Các tác giả đã thực hiện việc thay các số hạng phi tuyến bằng



các số hạng phi tuyến bậc cao hơn trước khi tuyến tính hóa với qui trình như sau





   kx

 g x





  k g x 1

 k g x 2

 g x x

k k ,

,

k  được tối ưu theo tiêu chuẩn cực tiểu hóa thông thường

(72)

2

2

2

g x ( )





k 1

min k 1

g x ( ) x

  

  

2

2

(73)

g x ( )





k 1

  k g x 2

min k

2

g x ( ) x

  

  

2



k x





 k g x 2





min k

Từ (73) ta có

3

3

/

x

/

x

 g x



,

,



k



k 1

2

k 1

4

2

2

/

x

  g x   g x

 g x



Các hệ số  1

65

3

3

/

x

/

x

k



4

2

2

  xg x 2 x

/

x

  g x   g x

  g x   g x

(74)

Tiêu chuẩn (73) còn được gọi là tiêu chuẩn điều chỉnh một bước. Áp dụng

tiêu chuẩn (73) cho hệ Duffing và Van der pol đã thu được sự cải thiện đáng kể

về độ chính xác của nghiệm. Elishakoff vcs (2009) [33] phát triển tiêu chuẩn của

2





Anh và Di Paola từ điều chỉnh một bước thành điều chỉnh hai bước với  qui trình











kx

  g x







 k g x 1

 k g x 2

  k g x 3

 k g x 4

 g x x

 g x x

  g x x

  

  

(75)

,

,

k  được xác định  theo  tiêu  chuẩn  cực  tiểu hóa

k k , 1

2

k k , 3

4

trong  đó các hệ  số

thông thường. Sau đó, các tác giả áp dụng qui trình (75) cho dao động Lutes and

Sarkani thu được kết quả tốt hơn so với qui trình điều chỉnh một bước ở trường

hợp mức độ phi tuyến thay đổi khá lớn.

3.2 Tiêu chuẩn đối ngẫu có trọng số

Tiêu chuẩn cực tiểu hóa bình phương trung bình của sai số phương trình do

Caughey  đề  xuất  thường  được  coi  là  cơ  sở  để  so  sánh  cho  các  phát  triển  lý

thuyết  của  phương  pháp  tuyến  tính  hóa  tương  đương.  Do  đó,  tiêu  chuẩn  của

Caughey còn được biết đến với các tên gọi khác  như  “tiêu chuẩn (hay phương

pháp)  kinh  điển”  hoặc “tiêu chuẩn  (hay phương  pháp)  thông  thường”. Các  kết

quả nghiên cứu trong (Roberts vcs 1990 [74], Socha 2008 [77], Proppe vcs 2003

điển cho kết quả tốt (sai số chấp nhận được) khi áp dụng cho các hệ có mức độ

phi tuyến yếu và độ chính xác giảm khi mức độ phi tuyến tăng lên.

Nguyễn Đông  Anh  (Anh  2010)  [17]  đề xuất  cách  tiếp  cận  mới  được  biết

với  tên gọi  cách  tiếp cận đối  ngẫu với  quan  điểm tạo  ra  một  sự  hài  hòa  trong

nghiên cứu, cho phép phát hiện bản chất của vấn đề một cách đầy đủ hơn. Tiêu

chuẩn đối ngẫu đầu tiên là tiêu chuẩn tuyến tính hóa có điều chỉnh do Anh và Di

Paola đề xuất hay còn được gọi là điều chỉnh một bước trong đó các tác giả đề

[70], Langley 1988 [49], Elishakoff vcs 2009 [33]) chỉ ra rằng tiêu chuẩn kinh

66

nghị  làm  tăng  mức  độ  phi  tuyến  của  hệ  ban  đầu  trước  khi  tuyến  tính  hóa.

Elishakoff  vcs  [2009]  [33]  phát  triển  tiêu  chuẩn  của  Anh  và  Di  Paola  từ  điều

chỉnh  một  bước  thành  điều  chỉnh  hai  bước  và  áp  dụng  cho  dao  động  Lutes

Sarkani có mức độ phi tuyến thay đổi mạnh. Tiêu chuẩn thứ hai là tiêu chuẩn sai

số bình phương trung bình toàn cục-khu vực được Anh vcs (2012c) [16] cải tiến

từ  tiêu chuẩn sai số bình phương trung bình khu  vực  trong đó  các tác giả phát

triển từ  miền đáp ứng tập trung sang miền đáp ứng toàn cục.  Anh (2010) [17]

nghiên cứu sự thay thế đối ngẫu từ hệ tuyến tính hóa trở về hệ phi tuyến ban đầu

đối với tiêu chuẩn kinh điển và áp dụng cho dao động ngẫu nhiên Duffing, sau

đó tác giả đã đề xuất một tiêu chuẩn đối ngẫu có trọng số tổng quát.

Ta có thể xem xét tiêu chuẩn đối ngẫu có trọng số được trình bày dưới dạng

2

2

như sau

S

  (1

p

)













dn

F nl

k F dn l

 p k F l

dn

 F nl

   min  , k

dn

(76)

p 

1/ 2

. Ta thấy rằng khi p trong đó p là trọng số được chọn trong khoảng  0

=0 sự thay thế lượt đi chiếm ưu thế hoàn toàn, khi p =1/2 thì sự thay thế lượt đi

và  thay  thế lượt về  (đối  ngẫu)  có vai  trò như  nhau.  Fnl là  hàm phi  tuyến,  Fl  là

dnk ,  là hệ số trở về. Các số hạng thứ

hàm tuyến tính với hệ số tuyến tính hóa

nhất và thứ hai của (76) tương ứng là các quá trình thay thế lượt đi và lượt về

được biểu diễn dưới dạng trung bình bình phương.

dnk  và hệ số trở về , tiêu chuẩn đối

ngẫu yêu cầu cực tiểu hóa tổng

dnS

dn

,

(77)



0

 S  k

dn

(78)



0

 dnS  

Từ (77), (78) ta có

,

(79)

k



dn

) 2

 p (1  pr 1

F F nl l 2 F l

Để xác định các hệ số tuyến tính hóa

67

2

2

trong đó ký hiệu

r



2

F F l nl 2 F F nl l

(80)

2r  trong công thức (80) mô tả quan hệ tương quan giữa hai hàm

Đại lượng

Fnl và Fl dưới dạng các mô men.

3.3. Những cải tiến của phương pháp đối ngẫu có trọng số

3.3.1. Cải tiến 1

Ta thấy rằng  hệ  số tuyến tính hóa tương  đương  trong (79) phụ thuộc vào

p 

1/ 2

. Việc xác định giá trị nào giá trị của trọng số p nằm  trong khoảng  0

của p là bài toán rất khó và là mục tiêu lâu dài của các nghiên cứu về tiêu chuẩn

đối ngẫu. Hiện nay,  hai giá trị giới hạn của p đang được sử dụng và quan  tâm

nhiều,  cụ thể  p=0 ứng với  tiêu  chuẩn  kinh  điển  và p  =1/2 đã được đề  nghị  và

đánh  giá  hiệu  quả  cho  các  hệ  phi  tuyến  ngẫu  nhiên  (Anh  vcs  2012a,  2012b)

[13,15].  Tuy  nhiên  các  áp  dụng  sơ  bộ  giá  trị  p=1/2  cho  các  hệ  phi  tuyến  dao

động tuần hoàn cho các kết quả không tốt, có thể là đối với các hệ dao động tiền

định cần có sự cân đối giữa các sự thay thế lượt đi và về. Qua nghiên cứu nhận

thấy  rằng  khi hệ  có  tính phi  tuyến  yếu  thì p xấp  xỉ 0,  tức là  tiêu  chuẩn  tương

đương kinh điển cho sai số nhỏ, và khi tính phi tuyến tăng thì p tăng theo. Trên

thực tế, với giá trị trọng số p=1/2 thì vai trò của sự thay thế lượt đi và lượt về là

thế thứ nhất trong khi sự thay thế đối ngẫu chỉ là thay thế thứ hai. Đối với hệ phi

tuyến vừa hai sự thay thế này không nên được cân bằng với nhau mà sự thay thế

thứ nhất nên có vai trò quan trọng hơn. Do đó luận án đề xuất kết hợp 2 giá trị

trọng số ở  biên bằng phép trung bình cộng. Cụ thể ta lấy 2 giá  trị của p trong

(79)

-  với p=0 ta có hệ số tương đương kinh điển (biên bên trái)

(81)

k



  0

dn

F F nl l 2 F l

tương đương. Tuy nhiên, ta cũng chú ý rằng việc thay thế kinh điển là sự thay

68

- với p=1/2 ta có hệ số tương đương đối ngẫu (biên bên phải)

k



dn

2

1 

2

r

F F nl l 2 F l

(82)

Với dạng  cải  tiến  thứ  nhất,  hệ  số  tương  đương  được  đề nghị  sẽ  là  giá  trị

2

trung bình cộng của hai hệ số tương đương biên trên

k



1





dn v _ 1

2

2

1 2

1 

2

r

  

  

F F nl l 2 F l

F F nl l 2 F l

 3 r   2 2 r



(83)

3.3.2. Cải tiến 2

Cũng  với lý luận như trên, ta cần sử dụng một phép trung bình nào đó với

trọng số ở 2 biên. Trong dạng cải tiến này, ta sử dụng trung bình tích phân chứ

không phải trung bình đại số như dạng cải tiến 1. Dạng cải tiến này cũng có thể

được lý luận theo các thuật ngữ toàn cục – địa phương được sử dụng trong (Anh

vcs 2013, 2014) [14] [65]. Theo cách lý luận này, ta có thể xem xét rằng các giá

trị của trọng số thay đổi trên một miền toàn cục của việc lấy tích phân. Và hệ số

tuyến tính hóa tương đương được đề xuất là giá trị trung bình của tất cả các hệ

số địa phương theo nghĩa tích phân. Cụ thể, ta lấy tích phân (79), hệ số được đề

1/2

nghị bằng:

k



k



1/ 2 2

dp



 p dp

dn v

_ 2

dn

) 2





1 1 / 2

 p (1  pr 1

0

0

F F nl l 2 F l

(84)

2

r

Tích phân này dẫn tới

 2 1

2

(85)

k





dn v

_ 2

 4

1 2 r

r

r 2

F F nl l 2 F l

 ln 1  

  

   

   

3.3.3. Cải tiến 3

Ta chú ý rằng trong tiêu chuẩn đối ngẫu (76) chứa 2 thành phần thể hiện sai

số giữa hàm phi tuyến Fnl và hàm tuyến tính Fl trong quá trình lượt đi và lượt về.

Tuy nhiên sau 2 lượt đi và về thì hàm phi tuyến sẽ bị thay đổi một hệ số  và sẽ

có sai số do  khác 1. Như vậy cách sử dụng 2 thành phần sai số chưa thực sự

69

đầy đủ.  Trong  dạng  cải  tiến  này,  3 dạng  sai  số  sẽ  được  xem  xét  với  trọng  số

2

2

2

bằng nhau. Cụ thể tiêu chuẩn đối ngẫu được bổ sung để trở thành:

S

























dn

F nl

k F dn l

k F dn l

 F nl

F nl

 F nl

   min  , k

dn

(86)

dn

Từ điều kiện cực tiểu



0,



0

 S  k

 S dn  

dn

(87)

Ta thu được biểu thức của hệ số tuyến tính hóa tương đương dạng cải tiến

3.

k



dn v

_ 3

2

3 

r

4

F F nl l 2 F l

(88)

3.4 Áp dụng cho dao động tự do của hệ phi tuyến dạng Duffing bậc cao

 x



x  2 1 n



0

Ta xét dao động tự do trong hệ phi tuyến dạng Duffing như sau:

(89)

với điều kiện đầu cho dịch chuyển

x

   (0) 1,

(90)

và cho vận tốc

 x

 (0) 0

(91)

Đây là hệ bảo toàn do đó ta có tần số dao động chính xác được tính theo

1

1

dx

n



(92)



,



4



4

 e

T e

1

dx 1  2 2 n





2

n

 2 T e

1



x

0

0

[2

x

 1 dx

1/ 2 ]



x

Phương trình tuyến tính tương đương với (89) sẽ có dạng

(93)

 x

kx



0

Áp dụng các công thức (79), (83), (85), (88) ta sẽ có các biểu thức tính tần

số riêng trong các trường hợp như sau:

- Tính theo tiêu chuẩn kinh điển (p=0)

công thức sau (Rao 2010)

70

2

n



2

2

n



2

x

t ( )

cos

)

t



k







2  kd

kd

2 x t ( )

2 cos (

 ( kd t )

F F nl l 2 F l

 kd

(94)

2

n



2

2

n



2

x

t ( )

cos

)

t



k







,

2  dn

dn

2

2

2

1 

(2

r

)

1 

(2

r

)

(2

r

)

1 

2 x t ( )

2 cos (

 ( dn t )

F F nl l 2 F l

 dn

- Tính theo tiêu chuẩn đối ngẫu trọng số với p=1/2

2

2

n



2

cos

t

)

2

r



4

n



2

 ( dn 2 ) cos (

cos

t

t

)

 ( dn

 dn

(95)

2

n



2

2

2

x

t ( )



k





2  dn v _ 1

dn v _ 1

2

2

2 x t ( )

F F nl l 2 F l

r  3   2 2 r



r  3   2 2 r



2

n



2

2

cos

)

 (

- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 1:



,

2

2 cos (



t _ 1 dn v t ) _ 1 dn v

r  3   2 2 r



2

2

n



2

cos

 (

)

2

r



4

n



2

t _ 1 dn v 2 cos (

cos

 (

)



)

t _ 1 dn v

t _ 1 dn v

(96)

- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 2:

71

2

2

r

 2 1





k









2  dn v

_ 2

dn v

_ 2

 4

1 2 r

r

r 2

F F nl l 2 F l

 ln 1  

  

   

   

2

2

n



2

2

r

x

t ( )

 2 1









 4

1 2 r

r

r 2

2 x t ( )

 ln 1  

  

   

   

2

2

n



2

2

r

cos

 (

t

)

 2 1





,

 4

1 2 r

r

r 2

2 cos (



_ 2 dn v t )

dn v

_ 2

  ln 1 

  

   

   

2

2

n



2

cos

 (

t

)

2

r



4

n



2

_ 2 dn v 2 cos (

cos

 (

t

)

)

dn v

_ 2

dn v t

_ 2

(97)

2

n



2

2

n



2

x

t ( )

cos

 (

)



k







,

2  dn v

_ 3

dn v

_ 3

2

2

2

3 

r

4

3 

r

4

3 

r

4

2 x t ( )

t _ 3 dn v t )

2 cos (



F F nl l 2 F l

dn v

_ 3

- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 3:

2

2

n



2

cos

 (

t

)

2

r



4

n



2

dn v _ 3 2 ) cos (

cos

 (



)

t _ 3

dn v

t _ 3

dn v

(98)

2 /

 

dt

    2

0

Sử dụng phép trung bình trong trường hợp dao động điều hòa:

(98). Đoạn mã MATLAB để tính các biểu thức được cho trong phụ lục.

So sánh các tần số tính toán với tần số chính xác được thể hiện trên Bảng 4.

trong đó  là tần số dao động ta sẽ tính được các biểu thức (94), (95), (96), (97),

Bảng 4. So sánh các tần số tính toán với tần số chính xác  dn v    n

dn v

dn v

_ 2

_ 3

_ 1

dn

kd

e

1  0.847  0.866  0.826  0.846  0.850  0.852

2  0.747  0.791  0.720  0.756  0.762  0.765

3  0.675  0.740  0.652  0.697  0.704  0.707

72

Kết quả trong bảng 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  cải

tiến cho kết quả tốt hơn tiêu chuẩn đối ngẫu 2 thành phần với các giá trị của n

không quá lớn (tính phi tuyến không quá lớn). Điều đó cho thấy tiêu chuẩn đối

ngẫu cải tiến có thể áp dụng tốt cho các hệ có tính phi tuyến trung bình.

3.5. Áp dụng cho dao động ngẫu nhiên

3

Chuyển sang dao động ngẫu nhiên ta xét hệ Duffing bậc ba như sau

 x



2

 hx



     x



x

  t

2 o

(99)

,

,

oh   là các hệ số dương và hệ số phi tuyến  sẽ thay đổi từ phi  ,

trong đó

  t

tuyến nhỏ sang phi tuyến trung bình, là quá trình ồn trắng Gauss với





0,

     t



  t

  t





  t

(100)



2

2

4

x

exp

x



 x

dx



2  o



1 2

1 4

2

  

  



Phương trình phi tuyến (99) có dịch chuyển trung bình bình phương là

x





c x

2

4

exp



x



 x

dx

2  o



h 4 2  4 h 2 

   1 2

1 4

  

  

  

     



(101)

Xét phương trình tuyến tính hóa tương đương ứng với (99)

 x



2

 hx



x



kx





  t

2  o

(102)

trong  đó  k là  hệ số  tuyến  tính hóa.  Theo  lý  thuyết  dao động  ngẫu  nhiên  tuyến

2

(103)

x





k

4

h

2   2  o



3

Áp dụng tiêu chuẩn đối ngẫu với

,

x , ta có

x

lF

nlF

2

2



x

F l

2

4

2

(104)





x



 3

x

F F nl l

3

6

2



2 

x



15

2 

x

2 F nl

tính, ta xác định theo công thức

73

2

2

2



3

x

2

Và

r





3

2

2

2

3 5

   1 5 

x

   x

(105)

Ta có hệ số tuyến tính hóa k dịch chuyển trung bình bình phương theo các

tiêu chuẩn như sau:

k





 3

2 x t ( )

kd

F F nl l 2 F l

- Tính theo tiêu chuẩn kinh điển

2

2

2

2

x





 3

x



x



2  o

kd

kd

kd

2  4 h

4

h



k

kd

2   2  o



(106)

k







2 x t ( )

dn

2

1 

2

r

15 7

F F nl l 2 F l

- Tính theo tiêu chuẩn đối ngẫu trọng số với p=1/2

2

2

2

2

x







x



x



2  o

dn

dn

dn

15 7

2  4 h

4

h



k

dn

2   2  o



(107)

2

k







2 x t ( )

dn v _ 1

2

18 7

F F nl l 2 F l

r  3   2 2 r

- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 1:

2

2

2

2

x







x



x



2  o

dn v _ 1

dn v _ 1

dn v _ 1

18 7

2  4 h

4

h

k

dn v _ 1



 2   2   o

2

2

r

(108)

 2 1



k







5



ln

2 x t ( )

dn v

_ 2

 4

r

r 2

1 2 r

20 3

7 10

  

  

  

   

F F nl l 2 F l

 ln 1  

  

   

   

(109)

2

2

2

2

x



ln

x



x



2  o

dn v

_ 2

dn v

_ 2

dn v

_ 2

20 3

7 10

2  4 h

4

h

k

  

  

   5  

   

dn v

_ 2

2   2   o



- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 3:

- Tính theo tiêu chuẩn đối ngẫu cải tiến dạng 2:  

74

k







2 x t ( )

dn v

_ 3

2

3 

r

4

45 17

F F nl l 2 F l

2

2

2

2

x







x



x



2  o

dn v

_ 3

dn v

_ 3

dn v

_ 3

45 17

2  4 h

4

h

k

dn v

_ 3

2   2   o



(110)

So  sánh  các  dịch  chuyển bình  phương  trung  bình  gần  đúng  và  chính xác

 h

1,



0.5,



1

   và hệ  số  phi tuyến   thay

o

được thể  hiện  trên  Bảng  5  với

2

2

2

2

2

x

x

đổi.

Bảng 5. So sánh các dịch chuyển bình phương trung bình  2 x x

x

x

dn v _ 1

dn v

_ 2

dn v

_ 3

cx

kd

dn



0,1  0.444  0.442  0.456  0.448 0.448 0.447

0,5  0.345  0.333  0.361  0.346 0.344 0.344

5 0.168  0.152  0.174  0.162 0.161 0.160

10  0.127  0.114  0.131  0.121 0.120 0.120

Các kết quả cho thấy xu hướng giống như trong ví dụ ở mục 3.4. Các dạng

cải tiến của tiêu chuẩn đối ngẫu đều cho kết quả tốt hơn tiêu chuẩn kinh điển.

Ngoài ra, kết quả sẽ tốt ở vùng có tính phi tuyến không quá nhỏ và quá lớn.

Kết luận chương 3

Trong chương này đã  trình bày tiêu chuẩn đối ngẫu có trọng số cho vấn đề

bằng không sẽ thu được tiêu chuẩn tuyến tính hóa kinh điển. Để giải quyết vấn

đề  chọn  giá  trị  trọng  số  như  thế  nào  đã  nghiên  cứu  đề  xuất  3  cách  lựa  chọn

tương ứng với 3 cải tiến. Áp dụng cho các hệ phi tuyến dạng đa thức là dạng hay

gặp trong bài toán ổn định flutter cho thấy trong miền khảo sát của các tham số

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 cho kết quả tốt hơn

tiêu chuẩn đối ngẫu 2 thành phần với các giá trị của n không quá lớn (tính phi

tuyến không quá lớn). Đ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 mạnh trung bình.

tuyến  tính  hóa  tương đương  hệ dao động phi tuyến, trong đó  khi  cho  trọng số

75

CHƯƠNG 4. ÁP DỤNG KỸ THUẬT TUYẾN TÍNH HÓA ĐỐI NGẪU

CHO BÀI TOÁN PHÂN TÍCH ĐÁP ỨNG PHI TUYẾN CỦA THIẾT

DIỆN CÁNH

4.1. Mô hình thiết diện cánh

Trong chương này ta sẽ khảo sát mộ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 trình bày trong Strganac vcs

2000. Mô hình thí nghiệm này được phát triển để cung cấp các số liệu đo trực

tiếp của đáp ứng khí động phi tuyến cũng  như  sự  tùy biến trong việc thay đổi

các điều kiện và tham số thí nghiệm. Độ lêch tâm của tâm khí động, khối lượng

của các thành phần khác nhau của hệ thống, mô men tĩnh, mô men quán tính của

thiết diện cánh, các đặc tính độ cứng và hình học của thiết diện cánh đều có thể

dễ dàng thay đổi trong các nghiên cứu khảo sát. Ngoài ra mô hình còn có phần

đuôi  cánh  có  thể  điều  khiển  để  kiểm  tra  các  thuật  toán  điều  khiển  khác  nhau

Trục đàn hồi

Dây cung trung bình

U

(hình 36).

Hình 36. Mô hình khí động với các bậc tự do tịnh tiến và xoắn

76

Hình 37. Mô hình giản đồ của bố trí thí nghiệm

Thiết diện cánh được gắn trên bệ đỡ linh động cho phép chuyển động với 2

bậc tự do. Thiết diện cánh còn bao gồm phần đuôi trải dài trên toàn bộ sải cánh

Đo quang uốn xoắn

Đo gia tốc uốn xoắn

Điều khiển  động cơ

Chíp LS7166

Điều kiện số

D/A

I/O số

Dịch dữ liệu

Chuyển đổi áp suất

A/D

Máy tính I, 90 MHz

Máy tính II, 333 MHz

NI PCI-MIO-16XE-10

(Hình 37, 38).

Hình 38. Mô hình giản đồ của số liệu đầu vào và đầu ra cho việc đo

và điều khiển mô hình thí nghiệm

77

Chuyển động tịnh tiến (plunge) được cung cấp bởi giá đỡ tịnh tiến. Chuyển

động  xoắn (pitch) độc lập  với chuyển động tịnh tiến được cung cấp bánh cam

gắn với giá đỡ. Thông thường thì chỉ có khối lượng của cánh cần phải phân tích

nhưng trong mô hình thí nghiệm thì do giá đỡ cũng chuyển động cùng với cánh

nên khối lượng tổng cộng của toàn bộ cánh và kết cấu đỡ phải được quan tâm

trong các phân tích. Thêm vào đó, mô men quán tính gắn với chuyển động xoắn

thì  không  bao gồm  đóng  góp  của  giá  đỡ  bởi  vì  giá  đỡ không có  chuyển động

xoắn.

Mô hình hệ thống đỡ cho phép tùy biến trong các điều kiện thí nghiệm và

các tham số. Độ cứng kết cấu được tạo ra bởi cặp bánh cam được thiết kế cho

phép thay đổi độ cứng tuyến tính và phi tuyến. Dạng của bánh cam, độ cứng và

độ kéo dãn ban đầu của lò xo tạo ra tính phi tuyến. Đáp ứng của hệ thống được

đo bằng các cảm biến gia tốc và các đầu đo quang. Vận tốc dòng khí được xác

định từ ống thổi gió.

Các ký hiệu của mô hình được mô tả trên Hình 36, trong đó nửa chiều dài

dây cung được ký hiệu bằng b, trục đàn hồi của mô hình có vị trí ở khoảng cách

ab so với điểm giữa dây cung. Khối tâm có vị trí ở khoảng cách xαb tính từ tâm

đàn hồi. Chuyển động tịnh tiến h là chuyển động thẳng đừng tại tâm đàn hồi và

có chiều dương hướng xuống. Góc xoắn a quanh trục đàn hồi có chiều dương

khi mũi cánh hướng lên. Góc xoay  của đuôi điều khiển có chiều dương hướng

xuống. Phương trình chuyển động của mô hình trên hình 36 đã được viết trong

vcs 2004 [68], Ko vcs 2002 [47], A. Abdelkefi vcs 2012 [10]) và có dạng sau:







 a



 









a





m h m x b a

W

T

 c h k h h

h

2 bsU c a l

2  U bsc  l

 h U

1 2

 a b U

  

  

 a  

  



(111)









  a a

 a I a

 a m x bh c a a

W

k a

2









a





2 b sU c m

a

2 2 U b sc m

 

 h U

1 2

 a b U

  

  

 a  

  

rất nhiều tài liệu tham khảo (Strganac vcs 2000 [78], Li vcs 2011 [53], Platanitis

78

trong đó mW là khối lượng của thiết diện cánh, mT là tổng khối lượng của cánh

cộng với giá đỡ, Ia là mô men quán tính quanh trục đàn hồi, kh là độ cứng tuyến

tính  của  chuyển  động  tịnh  tiến  trong  khi  ka  là  độ  cứng  phi  tuyến  của  chuyển

động xoắn, ch và ca là các hệ số cản, U là vận tốc dòng,  là khối lượng riêng

không khí, s là sải cánh, cla và cl là các hệ số khí động nâng trên một đơn vị góc

xoắn,  cma  và cm  là  các  hệ số  khí  động mô men  trên một đơn  vị  góc  xoắn.  So

sánh phương trình (111) với (24) ta chỉ ra một số thay đổi sau:

- Ký hiệu ah được chuyển thành a để phù hợp với dạng của các phương trình

trong các tài liệu tham khảo đã nêu.

- Mô men tĩnh S được thay cụ thể bằng mWxab, bằng khối lượng cánh nhân với

khoảng cách từ khối tâm đến tâm đàn hồi.

- Có sự xuất hiện thêm các hệ số cản ch và ca.

- Có sự xuất hiện thêm độ dài sải cánh s vì ở đây ta tính trên toàn bộ cánh chứ

không phải một dải có độ dài đơn vị.

- Xuất hiện thêm thành phần lực và mô men do ảnh hưởng của góc điều khiển

 ở đuôi cánh.

- Độ cứng xoắn ka bây giờ  không còn là hằng  số mà là một hàm  của a,  thể

hiện tính phi tuyến đối với bậc tự do xoắn, là hiện tượng phi tuyến hay gặp nhất

trong bài toán thiết diện cánh.

đa thức bậc 4 (Strganac vcs 2000):

(112)











k a

k a 1

2 a a k k 2 a a 3

3 a k 4 a

4 a k 5 a

Chú ý rằng nếu ka là hằng số thì ta hệ tuyến tính.

Trong mô hình được trình bày, độ cứng ka được xét là một hàm của a dưới dạng

4.2. Phương trình xác định vận tốc tới hạn

Như đã trình bày trong chương 1, vấn đề điều khiển thiết diện cánh là bài

toán lớn và được nhiều nhà nghiên cứu quan tâm với các thuật toán điều khiển

đề xuất khác nhau. Trong luận án này, do tập trung chính vào cách tiếp cận đối

ngẫu nên sẽ chỉ quan tâm tới thuật toán điều khiển đơn giản là điều khiển PID.  79

Đây  là  thuật  toán  điều  khiển  đơn  giản  nhất  và  được  sử  dụng  nhiều  nhất.  Góc

điều khiển sẽ được xác định là tổng của thành phần tỷ lệ (P:proportional), thành

phần  tích  phân  (I:Integral)  và  thành  phần  vi  phân  (D:derivative)  nhân  với  một

vài hệ số khuếch đại nào đó. Như vậy bộ điều khiển PID được viết dạng:

 a K

 dt K

K





   a d

p

i

 a trong  đó  Kp,  Ki  và  Kd  là  các  hệ  số  khuếch  đại.  Thay  (113)  vào  phương  trình

(113)

(111) ta có thể đưa phương trình về dạng phương trình trạng thái giống như (41)

1 0 0 0 1 0 0 0 1

0 0 0

0 0 0



0 0 0

0 0 1

0 0 0

1 0 0

0 1 0

1

2

3

4

0 0 0 0 0 0

k 1 k

k k

6

m m 2 m m 3

2

7

k k 8

k k 9

k 5 k 10

       

       

       

       

   h    a     a    h      a  

       

h   a    a    h   a 

(114)

trong  đó các hệ  số mi  (i=1,..3),  ki  (i=1,..10)  được  xác định  theo  các  biểu  thức

 m m m m x b m I





;

;

W

a

3

T

2

1

; a

2

 

 





;

 



U bsc K

;

k 1

k k ; h

2

3

 l

i

 2 bsU c a l

 c K k  l p

k



k

 





 a b c UK

;

4

   c h

bsUc ; a l

5

bsU c a l

 l

d

1 2

  

  

  

  

như sau:

2

2

k



0;

k

 







;

6

7

k a

a

c K  m

p



2 2

2





;





k 8

U b sc K k m



i

 b sU c m b sUc m

9

; a

2

 







k 10

c a

b sU c m

a

 a b Uc m



dK

1 2

  

  

  

  

Sau khi xác định được các hệ số thì bài toán flutter có thể giải được bằng

các phương trình (44), (45). Sử dụng các biểu thức (115) vào (45) ta sẽ thu được

phương trình một biến đối với U được viết lại (để dễ theo dõi) như sau:

(115)

80

2







k k 1 8

k k 6 3

1

2 2

















 k m m k m k 4

2 5

3

1 10

k k 6 2

k k 3 9

k k 8 4

k k 1 7

m k 2 9

   

   



















k k 1 7

k k 3 9

k k 8 4

k k 7 4

k k 1 10

k k 2 9





  m m m 3     















k k 6 2 

k k m k m k 6 5 1 8 2 3 



k k 4 10

2 2

2 6

k k 5 9

k m 1 3

k k 6 3

k k 1 8

   















k k 1 10

2 2

3

k k m k m k 6 5 1 8 2

3

k k 2 9

















     m m m k k 1 7 4   

 



 m k m k m k 1 7  k m m k m k 4

m k 2 9

2 5

1 10

3

k k 4 10

m k m k m k 1 7

2 2

2 6

k k 5 9

k m 1 3

   

   

(116)

Phương trình này có thể dễ dàng giải bằng MATLAB nhờ lệnh fzero. Sau

khi tính được vận tốc tới hạn thì tần số dao động LCO được xác định từ (44) (và









k k 4 10

2 2

k k 5 9

k m 1 3





 



1

1 10







2 5 

k k 7 4

k k 2 9

m k 2 9 k k 6 2

k k 1 7





 

 

 

 

 m k m k m k 2 6 1 7 2  m m m 2 3 k k m k m k 6 5 1 8 2 3 k m m k m k 4

k k 1 10 m k 2 9

2 5

3

1 10

 k k k k 6 3 1 8 k m m k m k 3 4  k k k k 3 9 8 4 2  m m m 2 3

1

   2    

     

     

     

được viết lại ở đây để tiện theo dõi):

(117)

Đoạn mã mô tả phương trình được thể hiện trong phụ lục.

Ta chú ý rằng các biểu thức được lập ở đây là đối với bài toán tuyến tính

(ka là hằng số). Trong bài toán phi tuyến thì ta sẽ dùng các phương pháp tuyến

tính hóa thích hợp để tính độ cứng tương đương của ka.

4.3. Áp dụng kỹ thuật tuyến tính hóa đối ngẫu

được tóm tắt ở đây trước khi áp dụng vào bài toán thiết diện cánh cụ thể. Gọi Fnl

là hàm phi tuyến cần được tuyến tính hóa và Fl là hàm tuyến tính, ký hiệu <•>

thể hiện phép lấy trung bình (trong chương này sẽ sử dụng phép trung bình điều

hòa). Xét ký hiệu:

2

2

(118)

r



2

F F l nl 2 F F nl l

Các  kết  quả  của  các  kỹ  thuật  đối ngẫu  được  trình  bày  trong  chương 3 sẽ

81

Hệ  số  tuyến  tính  hóa  tương  đương  được  tính  theo  từng  trường  hợp  như

sau:

* Trường hợp tuyến tính hóa kinh điển

k



kd

F F nl l 2 F l

,                                               (119)

* Trường hợp tuyến tính hóa đối ngẫu thông thường

k



dn

2

1 

2

r

F F nl l 2 F l

,                                        (120)

2

* Trường hợp tuyến tính hóa đối ngẫu cải tiến 1

k



dn v _ 1

2

F F nl l 2 F l

 3 r   2 2 r



(121)

2

r

 2 1

2

* Trường hợp tuyến tính hóa đối ngẫu cải tiến 2

k





dn v

_ 2

 4

1 2 r

r

r 2

F F nl l 2 F l

 ln 1  

  

   

   

(122)

* Trường hợp tuyến tính hóa đối ngẫu cải tiến 3

k



dn v

_ 3

2

3 

r

4

F F nl l 2 F l

(123)

Trong trường hợp độ cứng phi tuyến như  trong công thức (112),  độ cứng

phi tuyến sẽ được tuyến tính hóa dưới dạng:

k a

  k a 1

k a

td

trong đó katd là độ cứng tương đương tính theo một trong các công thức (119),

(120), (121), (122), (123). Đối với độ cứng phi tuyến như trong công thức (112)

thì hàm phi tuyến và hàm tuyến tính là:









   a

F nl

2 a a k k 2 a a 3

3 a k 4 a

F l



 4 a a k ; 5 a

Giả sử đáp ứng điều hòa của góc xoắn có dạng:

a



sinA

 t

(124)

82

với  là tần số dao động. Sử dụng phép trung bình trong trường hợp dao động

2 /

 

dt

    2

0

điều hòa:

2

5

3



2

4 A k a 5 8

2

  

  

r





4

2

2

5

3

F F l nl 2 F F nl l

+

+

+

k k a a 3 5

2 k a 4

k k a a 2 4

2 k a 3

2 k a 5

2 k a 2

2 A k a 3 4 35 

 + 2

 + 2



6 A 64

8 A 63 128

A 8

A 4

3

5





2 A k a 3 4

4 A k a 5 8

F F nl l 2 F l

ta sẽ tính được các biểu thức:

Vậy  hệ  số  tuyến  tính  hóa  tương  đương  (124)  trong  các trường  hợp  sẽ  có

dạng:

3

5

* Trường hợp tuyến tính hóa kinh điển





k a

kd

a k

1

2 A k a 3 4

4 A k a 5 8

,                                   (125)

3

5

* Trường hợp tuyến tính hóa đối ngẫu thông thường







k a

dn

k a 1

2

1 

2

r

2 A k a 3 4

4 A k a 5 8

  

  

, (126)

2

3

5

* Trường hợp tuyến tính hóa đối ngẫu cải tiến 1







k a

dn v _ 1

k a 1

2

2 A k a 3 4

4 A k a 5 8

  

  

 3 r   2 2 r



* Trường hợp tuyến tính hóa đối ngẫu cải tiến 2

2

r

3

5

(127)

 2 1

2

(128)











k a

dn v

_ 2

k a 1

 4

1 2 r

r

r 2

2 A k a 3 4

4 A k a 5 8

 ln 1  

  

   

         

* Trường hợp tuyến tính hóa đối ngẫu cải tiến 3

3

5

(129)







k a

dn v

_ 3

k a 1

2

3 

r

2 A k a 3 4

4

4 A k a 5 8

  

  

83

Thay các hệ  số tuyến  tính  hóa  tương  đương (125)  hoặc  (126)  hoặc  (127)

hoặc (128) hoặc (129) vào vị trí của ka trong (115), rồi thay vào (116) ta sẽ thấy

vận tốc U là hàm của biên độ dao động xoắn A. Từ đó ta sẽ vẽ được đường đặc

trưng biên độ LCO với vận tốc tới hạn. Biết được quan hệ giữa U và A thì biểu

thức (117) cũng cho ta mối quan hệ giữa tần số LCO với vận tốc tới hạn.

4.4. Các ví dụ và tính toán bằng phương trình vi phân

Trong mục này ta sẽ tóm tắt các số liệu của các trường hợp ví dụ sẽ được

nghiên cứu cũng như các kết quả tính toán số cho các trường hợp.

4.4.1. Số liệu đầu vào

Các số liệu của các ví dụ lấy theo tài liệu Li vcs 2011 [53]. Ta sẽ khảo sát 5

ví dụ sau đây.

Ví dụ 1: Thiết diện cánh phi tuyến bậc 5, không có điều khiển.

Các số liệu của các tham số trình bày trong mục 4.1 được cho trên bảng 6.

a b (m)

Bảng 6. Số liệu của ví dụ 1  mT (kg)

mW (kg) Ia (kg.m2)

-0.6847 0.135 12.387 2.049 0.0558

kh (N/m) ch (kg/s) xa cla cma

0.3314 2884.4 27.43 6.38 -1.16

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

0

0

s (m)   (kg/m3)  ca (kg.m2/s)  ka1 (N.m)  ka2 (N.m)

Các số liệu được cho trên bảng 7.

Ví dụ 2: Thiết diện cánh phi tuyến bậc 3, không có điều khiển.

84

a b (m)

Bảng 7: Số liệu của ví dụ 2  mT (kg)

mW (kg) Ia (kg.m2)

-0.6719 0.1905 15.57 5.23 0.1419

kh (N/m) ch (kg/s) xa cla cma

0.5721 2844 27.43 6.757 -1.162

s (m)   (kg/m3)  ca (kg.m2/s)  ka1 (N.m)  ka2 (N.m)

0.5945 1.225 0.0184 12.77 53.47

ka3 (N.m)  ka4 (N.m) ka5 (N.m) cl cm

1003 0 0 0 0

Ví dụ 3: Thiết diện cánh phi tuyến bậc 5, điều khiển tỷ lệ (P control).

Các số liệu được cho trên bảng 8.

Bảng 8: Số liệu của ví dụ 3

b (m) a mT (kg) mW (kg) Ia (kg.m2)

-0.6847 0.135 12.387 2.049 0.0558

kh (N/m) ch (kg/s) cla cma xa

0.3314 2884.4 27.43 6.38 -1.16

s (m)  (kg/m3)  ca (kg.m2/s)  ka1 (N.m)  ka2 (N.m)

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)

Thay đổi để khảo sát

0

0

1.225 0.036 6.833 9.967 0.6

Ví dụ 4: Thiết diện cánh phi tuyến bậc 5, điều khiển tích phân (I control)

Các số liệu được cho trên bảng 9.

85

Bảng 9: Số liệu của ví dụ 4

a b (m) mT (kg) mW (kg) Ia (kg.m2)

-0.6847 0.135 12.387 2.049 0.0558

kh (N/m) ch (kg/s) xa cla cma

0.3314 2884.4 27.43 6.38 -1.16

s (m)  (kg/m3) ca (kg.m2/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.

a b (m)

Bảng 10: Số liệu của ví dụ 5  mT (kg)

mW (kg) Ia (kg.m2)

-0.6847 0.135 12.387 2.049 0.0558

kh (N/m) ch (kg/s) xa cla cma

s (m)   (kg/m3)

ca (kg.m2/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

0.3314 2884.4 27.43 6.38 -1.16

86

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

h

  0 ,

h a    0 ,

  0

đầu khác đượ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

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.

tốc tới hạn có giá trị nhỏ nhất.

87

Để đả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

0.35

0.25

h(0)=0 h(0)=0.01 h(0)=0.05

0.15

0.05

là kết quả của luận án còn hình 40 là kết quả trong Li vcs 2011 [53].

d a r , ) 0 (

-0.05

a

-0.15

-0.25

-0.35

7.5

8

8.5

10.5

11

11.5

12

9

10

9.5 Vận tốc tới hạn , m/s

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)

88

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í

TTH kinh điển

) s / m

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

TTH đối ngẫu  cải tiến 3 TTH đối ngẫu  thông thường

(   n ạ h   i ớ t   c ố t   n ậ V

12.5 12 11.5 11 10.5 10 9.5 9 8.5 8 7.5

0.09

0.1

0.11

0.13

0.14

0.12 Biên độ LCO (rad)

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.

Hình 41: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 1  89

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 Vận  tốc  tới  hạn

(rad) (m/s)

Tính toán số 11.66 7.94

12.25 (sai số TTH kinh điển 7.9455 5.06%)

TTH đối ngẫu thông 10.5666  (sai  số 7.9455 9.38%) thường

0.1484 TTH  đối  ngẫu  cải 11.3431  (sai  số 7.9456 2.72%) tiến 1

TTH  đối  ngẫu  cải 11.51 (sai số 7.9457 1.29%) tiến 2

TTH  đối  ngẫu  cải 11.591 (sai  số 7.9456 0.59%) tiến 3

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ể

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.

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

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.

90

TTH kinh điển

) s / m

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

TTH đối ngẫu  cải tiến 3 TTH đối ngẫu  thông thường

(   n ạ h   i ớ t   c ố t   n ậ V

11.5 11.4 11.3 11.2 11.1 11 10.9 10.8 10.7 10.6 10.5

0.09

0.11

0.13

0.15

0.17

Biên độ LCO (rad)

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 Vận  tốc  tới  hạn

(rad) (m/s)

Tính toán số 10.525 11.29

11.4404  (sai  số TTH kinh điển 10.5246 1.33%)

thường

1.99%)

0.1743

TTH  đối  ngẫu  cải

11.2457  (sai  số

10.5246

tiến 1

0.39%)

TTH  đối  ngẫu  cải

11.2782  (sai  số

10.5247

tiến 2

0.1%)

TTH  đối  ngẫu  cải

11.2939  (sai  số

10.5246

tiến 3

0.03%)

TTH đối ngẫu thông 11.0654  (sai  số 10.5247

91

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à

11

TTH kinh điển

10.8

) s /

m

10.6

(

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

10.4

n ạ h

10.2

10

i ớ t   c ố t

TTH đối ngẫu  cải tiến 3 TTH đối ngẫu  thông thường

9.8

n ậ V

9.6

9.4

0.07

0.075

0.085

0.09

0.095

0.08 Biên độ LCO (rad)

bảng 15, 16.

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

Vận  tốc  tới  hạn

(rad)

(m/s)

92

Tính toán số 10.465 9.415

10.9671  (sai  số TTH kinh điển 9.4207 4.8%)

TTH đối ngẫu thông 9.8347 (sai  số 9.421 6.02%) thường

0.0972 TTH  đối  ngẫu  cải 10.1836  (sai  số 9.4209 2.69%) tiến 1

TTH  đối  ngẫu  cải 10.2784  (sai  số 9.4207 1.78%) tiến 2

TTH kinh điển

) s /

m

(

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

TTH đối ngẫu  cải tiến 3

n ạ h

i ớ t   c ố t

TTH đối ngẫu  thông thường

n ậ V

15.5 14.5 13.5 12.5 11.5 10.5 9.5 8.5 7.5 6.5

0.1

0.12

0.14

0.16

0.2

0.22

0.18 Biên độ LCO (rad)

TTH  đối  ngẫu  cải 10.3285  (sai  số 9.421 1.30%) tiến 3

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

Vận  tốc  tới  hạn

(rad)

(m/s)

93

Tính toán số 14.745 6.985

15.1676  (sai  số TTH kinh điển 6.9877 2.87%)

TTH đối ngẫu thông 13.9926  (sai  số 6.988 5.1%) thường

0.2255 TTH  đối  ngẫu  cải 14.5797  (sai  số 6.9877 1.12%) tiến 1

TTH  đối  ngẫu  cải 14.6995  (sai  số 6.9879 0.31%) tiến 2

TTH  đối  ngẫu  cải 14.7562  (sai  số 6.9877 0.08%) tiến 3

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

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.

độ 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

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

94

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

10

TTH kinh điển

9.5

) s /

m

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

9

TTH đối ngẫu  cải tiến 3

8.5

(   n ạ h   i ớ t   c ố t

8

TTH đối ngẫu  thông thường

n ậ V

7.5

7

0.1

0.11

0.12

0.14

0.15

0.13 Biên độ LCO (rad)

thể hiện trên hình 46 và bảng 18.

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 Vận  tốc  tới  hạn

Tính toán số

9.385

7.055

9.5412  (sai  số

TTH kinh điển

7.058

1.66%)

0.1595

TTH đối ngẫu thông

8.848

(sai

số

7.079

thường

5.72%)

TTH  đối  ngẫu  cải

9.1912  (sai  số

7.0577

tiến 1

2.06%)

(rad) (m/s)

95

TTH  đối  ngẫu  cải 9.2598  (sai  số 7.0578 tiến 2 1.33%)

) d a r (

a   g n ộ đ   o a d

c ó G

0.01 0.008 0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01

0

20

40

80

100

120

60 Thời gian (s)

TTH  đối  ngẫu  cải 9.2924  (sai  số 7.0577 tiến 3 0.99%)

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 TTH đối TTH đối TTH đối TTH đối

điển ngẫu thông ngẫu cải ngẫu cải ngẫu cải

Các giá trị

-2.7795

-2.7843

-2.7819

-2.7818

-2.7817

riêng của

+14.1352i

+14.1308i

+14.1330i

+14.1331i

+14.1331i

(37) ứng với

-2.7795

-2.7843

-2.7819

-2.7818

-2.7817

trường hợp

-14.1352i

-14.1308i

-14.1330i

-14.1331i

-14.1331i

biên độ rất

14.1336i

14.1300i

14.1318i

14.1319i

14.1319i

nhỏ =0.01

-14.1336i

-14.1300i

-14.1318i

-14.1319i

-14.1319i

rad

thường tiến 1 tiến 2 tiến 3

1.6001

1.6078

1.6039

1.6037

1.6037

96

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

hình 47, 48 và bảng 19,20.

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

97

2.9

TTH kinh điển

2.88

) s / m

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

2.86

TTH đối ngẫu  cải tiến 3

2.84

TTH đối ngẫu  thông thường

(   n ạ h   i ớ t   c ố t   n ậ V

2.82

2.8

0.07

0.075

0.08

0.085

0.09

Biên độ LCO (rad)

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 Vận  tốc  tới  hạn

(rad) (m/s)

Tính toán số 2.87 2.805

2.8857  (sai  số TTH kinh điển 2.8063 0.55%)

thường

1.01%)

0.0932

TTH  đối  ngẫu  cải

2.8606  (sai  số

2.8063

tiến 1

0.33%)

TTH  đối  ngẫu  cải

2.8649  (sai  số

2.8063

tiến 2

0.18%)

TTH  đối  ngẫu  cải

2.867

(sai

số

2.8063

tiến 3

0.1%)

TTH đối ngẫu thông 2.841 (sai số 2.8063

98

9.85

9.8

TTH kinh điển

) s / m

9.75

- Dấu X: kết quả số - Nét đứt: đường  cong không ổn định

9.7

TTH đối ngẫu  cải tiến 3

9.65

9.6

TTH đối ngẫu  thông thường

9.55

(   n ạ h   i ớ t   c ố t   n ậ V

9.5

9.45

0.05

0.055

0.065

0.06 Biên độ LCO (rad)

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 Vận  tốc  tới  hạn

(rad) (m/s)

Tính toán số 9.78 9.5

9.8426  (sai  số TTH kinh điển 9.4846 0.64%)

thường

2.24%)

0.0655

TTH  đối  ngẫu  cải

9.6472  (sai  số

9.4846

tiến 1

1.36%)

TTH  đối  ngẫu  cải

9.6703  (sai  số

9.4846

tiến 2

1.12%)

TTH  đối  ngẫu  cải

9.6826  (sai  số

9.4846

tiến 3

1%)

TTH đối ngẫu thông 9.5608  (sai  số 9.4846

99

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

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

Kết luận chương 4

100

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.

101

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

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

đã 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

102

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,

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

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ố

103

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.

104

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.

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.

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.

105

106

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

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ự.

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ĩ

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

107

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.

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.

16. Anh,  N.D.,  Hung,  L.X.,  Viet,  L.D.  (2012c)  Dual  approach  to  local  mean

108

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.

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.

27. Crandall  S.H.: (2006), "A  half-century  of stochastic  equivalent

109

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

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.

edition, Cambridge University Press.

110

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

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.

equivalent linearization method. J. Sound Vib.  127:271-281.

111

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–

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.

69.

112

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

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.

and  Monte-Carlo  simulation  in  stochastic  dynamics.  J.  Probab.  Eng.  Mech.

113

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

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.

alleviation using an electro-magnetic dry friction damper, part 2: experiment. J

114

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.

115

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:

2

m







  ...

p 0

p z 1

p z 2

p z m

2

n





  ...

  p z   q z

q 0

 q z q z 1 2

q z n

Cho p và q là hai đa thức khác không, có bậc tương ứng là m và n. Khi đó:

p

...

0 ... 0

Ma trận Sylvester tương ứng với p và q sẽ có kích cỡ (m+n)(m+n) và thu được



m

p m

 1

p 1

p 0

theo cách sau:  - Hàng đầu tiên là: 

-  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.

0 ... 0

...

q

- 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



q n

q 0

q 1

 1

n

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: 

- Các hàng tiếp theo thu được bằng cách như đã nói ở trên.

2







p 0

p z 1

p z 2

2





  p z   q z

q 0

 q z q z 1 2

Ma trận Sylvester tương ứng với 2 đa thức này là:

2

S



p q ,

0 p 0 0

p 2 0 q 2 0

q 0

p 1 p q 1 q 2

p 0 p 1 q 0 q 1

     

     

Xét ví dụ 2 đa thức bậc 2 có dạng:

116

2





















p qR

,

p q 2 0

q p 2 0

p q 0 1

p q 1 0

p q 2 1

q p 2 1

và kết thức của 2 đa thức có dạng:

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]

117

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;

118

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

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

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)/A1<1e-2 | A2>A1) & (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)

119