Nghiên cứu ảnh hưởng của độ cứng đến lực nâng của cánh vẫy kiểu côn trùng

Chia sẻ: Dạ Thiên Lăng | Ngày: | Loại File: PDF | Số trang:9

Thêm vào BST

Báo xấu

2
lượt xem 0
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Trong bài báo "Nghiên cứu ảnh hưởng của độ cứng đến lực nâng của cánh vẫy kiểu côn trùng", các tác giả xây dựng một chương trình tính toán khí động lực học đàn hồi cho cánh vẫy kiểu côn trùng dựa trên cách tiếp cận động lực học đa vật, trong đó kết cấu cánh được mô hình hóa dưới dạng hệ vật – lò xo và kết hợp với bộ giải khí động bằng phương pháp xoáy không dừng UVLM mở rộng. Mời các bạn cùng tham khảo!

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Nghiên cứu ảnh hưởng của độ cứng đến lực nâng của cánh vẫy kiểu côn trùng

212 Tuyển tập công trình Hội nghị Cơ học toàn quốc lần thứ XI, Hà Nội, 02-03/12/2022 Nghiên cứu ảnh hưởng của độ cứng đến lực nâng của cánh vẫy kiểu côn trùng Lê Vũ Đan Thanh1,*, Đặng Ngọc Thanh1 và Nguyễn Anh Tuấn1 1 Học viện Kỹ thuật Quân sự *Email: danthanh@lqdtu.edu.vn Tóm tắt. Cũng giống như các loài côn trùng trong tự nhiên, cánh của các thiết bị bay (TBB) cánh vẫy phỏng côn trùng thường được làm bằng các vật liệu nhẹ để giảm khối lượng. Khi vẫy trong không khí, chúng bị biến dạng và có ảnh hưởng lớn đến đặc tính khí động của TBB. Do đó, việc lựa chọn độ cứng phù hợp cho kết cấu cánh có ý nghĩa quan trọng trong quá trình thiết kế các TBB loại này. Trong bài báo này, các tác giả xây dựng một chương trình tính toán khí động lực học đàn hồi cho cánh vẫy kiểu côn trùng dựa trên cách tiếp cận động lực học đa vật, trong đó kết cấu cánh được mô hình hóa dưới dạng hệ vật – lò xo và kết hợp với bộ giải khí động bằng phương pháp xoáy không dừng UVLM mở rộng. Từ đó, xác định lực nâng trung bình trong mỗi chu kỳ với các độ cứng uốn khác nhau. Kết quả cho thấy, lực nâng lớn nhất đạt được khi tần số dao động riêng thứ nhất của kết cấu cánh xấp xỉ hai lần tần số vẫy. Từ khóa: Cánh vẫy kiểu côn trùng, Khí động lực học đàn hồi, Động lực học đa vật. 1. Mở đầu Các loài côn trùng có những đặc tính bay vượt trội như khả năng bay treo, bay lùi, cất – hạ cánh nhanh [1]. Cánh của chúng là kết cấu nhẹ, đàn hồi, thường là một tấm màng da được trợ lực bởi các gân hình côn, to ở gốc cánh và thon dần về mút cánh. Côn trùng hiện là nguồn cảm hứng cho nhiều nghiên cứu về TBB cánh vẫy phỏng sinh học. Cánh của TBB loại này phải làm bằng các vật liệu đủ nhẹ để đảm bảo về mặt khối lượng cất cánh, cũng như giảm công suất hoạt động của động cơ. Để tăng độ cứng của cánh, các kết cấu thường được gia cố thêm hệ thống khung carbon tương tự cấu trúc gân của cánh côn trùng [2]. Thực nghiệm cho thấy khi vẫy trong không khí, cánh bị biến dạng và có ảnh hưởng lớn đến các đặc tính khí động của TBB [3]. Các nghiên cứu về ảnh hưởng của độ cứng đến khả năng tạo lực của cánh vẫy kiểu côn trùng chủ yếu bằng thực nghiệm. Mountcastle và Combes [4] đã tiến hành thử nghiệm trên cánh của loài ong vò vẽ và thấy rằng những biến dạng thụ động của cánh có thể làm tăng lực nâng. Các nghiên cứu thực nghiệm của Campos và cộng sự [5] và Fu và cộng sự [6] lại cho thấy cánh quá mềm sẽ dẫn đến biến dạng uốn lớn, theo đó làm suy giảm các xoáy và giảm khả năng tạo lực. Cũng có một số nghiên cứu mô phỏng về vấn đề này, như nghiên cứu của Du và Sun [7] cho thấy sử dụng cánh mềm làm tăng lực nâng 10% và giảm công suất 5% so với cánh cứng. Truong [8] sử dụng chương trình tính toán tương tác kết cấu – chất lưu (FSI) để nghiên cứu 3 mô hình cánh vẫy, gồm cánh cứng, cánh mềm và cánh rất mềm. Kết quả cho thấy, cánh mềm có chất lượng khí động cao hơn so với cánh cứng, còn cánh rất mềm lại có chất lượng khí động kém hơn cánh cứng. Các tác giả cho rằng có một khoảng độ cứng của cánh mà tại đó tối ưu về mặt khí động. Tuy nhiên, do các mô hình mô phỏng FSI thường đòi hỏi khối lượng tính toán lớn, nên các nghiên cứu ở trên thường chỉ khảo sát được ít trường hợp, và kết quả cơ bản chỉ mang tính chất định tính. Trong bài báo này, các tác giả xây dựng một mô hình tương tác kết cấu – chất lưu FSI cho cánh vẫy kiểu côn trùng dựa trên cách tiếp cận động lực học đa vật. Kết cấu cánh được mô hình hóa dưới dạng hệ các vật cứng nối với nhau bằng lò xo uốn và xoắn. Hệ vật – lò xo được kết hợp với chương trình tính toán khí động lực học bằng xoáy không dừng UVLM mở rộng. Cách tiếp cận này cho phép giảm thời gian tính toán so với các phương pháp nghiên cứu FSI trước đây, do đó có thể áp dụng để khảo sát ảnh hưởng của độ cứng uốn đến lực nâng của cánh một cách cụ thể hơn.
213 Lê Vũ Đan Thanh, Đặng Ngọc Thanh, Nguyễn Anh Tuấn 2. Mô hình tính toán Bài báo này nghiên cứu mô hình cánh của loài nhậy Manduca Sexta. Phương pháp xây dựng mô hình kết cấu cánh đã được trình bày cụ thể trong một số nghiên cứu trước đây của các tác giả [9-11], sẽ được nêu lại một cách ngắn gọn để đảm bảo tính thống nhất của bài báo. Mô hình phần tử hữu hạn (PTHH) của cánh (Hình 1a) được xây dựng dựa trên dữ liệu thực nghiệm của O’Hara và Palazotto [12]. Với giả thiết rằng trục đàn hồi đi qua điểm mút cánh và gốc cánh, kết cấu cánh được mô hình hóa bằng một dầm Euler-Bernoulli chịu uốn và xoắn. Các đặc tính quán tính như khối lượng, mô men quán tính và vị trí trọng tâm được xác định từ mô hình PTHH cho các dải dọc theo sải cánh. Độ cứng chống uốn và chống xoắn của dầm được tính từ mô hình PTHH bằng cách sử dụng hai hệ điểm đo uốn và đo xoắn cách đều nhau (Hình 1b). Các điểm đo uốn nằm dọc theo trục đàn hồi, còn các điểm đo xoắn nằm gần mép sau cánh. Độ cứng chống uốn EI của dầm được xác định dựa trên độ cong đo tại mỗi điểm đo uốn khi tác dụng một lực điểm lên mút cánh theo hướng thẳng đứng. Độ cứng chống xoắn GJ được xác định dựa trên góc xoắn tại mỗi dải cánh khi lần lượt cho các lực điểm tác dụng lên các điểm đo xoắn [10]. Điểm đo xoắn Điểm đo uốn a) b) Hình 1. a) Mô hình phần tử hữu hạn của cánh loài nhậy Manduca Sexta; b) Hệ các điểm đo độ cứng Theo quả nghiên cứu thực nghiệm của Combes và Daniel [13], phân bố độ cứng dọc theo sải cánh có thể được xấp xỉ bằng hàm mũ có dạng: EI  a1b exp(  a2 r ) b (1) GJ  a1t exp(  a2 r ), t trong đó 𝑟̅ là tọa độ tương đối. Do cấu trúc phức tạp của cánh nên phương pháp xác định độ cứng như trên không thể đảm bảo chính xác hoàn toàn. Vì vậy, các hệ số độ cứng này cần được hiệu chỉnh qua quá trình tối ưu hóa. Hàm chi phí dùng cho tối ưu để làm giảm thiểu sự khác biệt giữa các dạng dao động riêng của mô hình dầm và mô hình PTHH. Ở đây chỉ xét đến hai dạng dao động riêng đầu tiên, là những dạng dao động có ảnh hưởng lớn nhất đến biến dạng cánh. Hàm chi phí tối ưu có dạng sau:  f1   f2  J( aib , ait )  1   1b  1  MAC1  1    2   2b  1  MAC 2  1  ,  f   f  (2)  FEM   FEM  b t i trong đó ai và ai tương ứng là hệ số thứ i của hàm xấp xỉ phân bố độ cứng uốn và xoắn; f FEM và fbi là tần số dao động riêng thứ i của mô hình PTHH và dầm; 1 và 2 là các trọng số. MAC1 và MAC2 là các số MAC (Modal Assurance Criterion) cho thấy sự tương đồng giữa các dạng dao động
214 Nghiên cứu ảnh hưởng của độ cứng uốn đến lực nâng của cánh vẫy kiểu côn trùng riêng của hai mô hình. Số MAC nhận giá trị từ 0 đến 1.0, trong đó số MAC bằng 1 tương ứng với trường hợp hai dạng dao động riêng hoàn toàn tương đồng [14]. Phân bố độ cứng và các hàm xấp xỉ tương ứng được biểu diễn trên hình 2. Các chấm tròn biểu diễn độ cứng đo được từ mô hình FEM, các đường nét đứt thể hiện hàm xấp xỉ độ cứng ban đầu, đường nét liền thể hiện phân bố độ cứng tối ưu. a) b) Hình 2. Phân bố độ cứng dọc theo sải cánh và hàm xấp xỉ tương ứng: a) độ cứng uốn; b) độ cứng xoắn Giá trị của các hệ số độ cứng được cho trong bảng 1. Bảng 1. Hệ số phân bố độ cứng của mô hình cánh 𝑎 𝑎 𝑎 𝑎 Trước khi tối ưu 8,65∙10-5 4,17 7,39∙10-5 5.77 Sau khi tối ưu 8,91∙10-5 4,43 7,56∙10-5 5,74 Bảng 2 thể hiện tần số và số MAC của các dạng dao động riêng của mô hình dầm trước và sau khi tối ưu, cũng như tần số dao động riêng của mô hình PTHH. Các dạng dao động riêng của mô hình PTHH và mô hình dầm được thể hiện trên Hình 3. Có thể thấy, tần số dao động riêng của mô hình dầm sau khi tối ưu và mô hình PTHH bằng nhau, số MAC giữa 2 mô hình cũng khá gần giá trị 1. Điều này chứng tỏ sự chính xác của mô hình cánh côn trùng dạng dầm. Bảng 2. Tần số và số MAC của các dạng dao động riêng f1 (Hz) f2 (Hz) MAC1 MAC2 J Mô hình PTHH 57,11 95,67 - - - Mô hình dầm trước khi tối ưu 58,88 94,72 0,96 0,81 0,12 Mô hình dầm sau khi tối ưu 57,11 95,67 0,95 0,80 0,10
215 Lê Vũ Đan Thanh, Đặng Ngọc Thanh, Nguyễn Anh Tuấn Hình 3. Các dạng dao động riêng của: a) mô hình PTHH, b) mô hình dầm Dầm phi tuyến với các thông số ở trên sẽ được xấp xỉ dưới dạng cơ hệ nhiều vật gồm N vật cứng có chiều dài bằng nhau được liên kết với nhau bằng các lò xo uốn và xoắn (Hình 4). Hệ được phát động bởi các chuyển động xoay của khớp cầu ở gốc cánh theo quy luật bay treo với tần số 26,1 Hz, được xác định từ thực nghiệm [15]. Vị trí và hướng tương đối của vật thứ k (k>1) so với vật thứ k- 1 được xác định bởi các góc uốn θk và xoắn αk của các lò xo nối các vật. Tọa độ suy rộng dùng để xác định trạng thái của hệ là Φ   2 ,  2 ,  3 ,  3 ,...,  N ,  N  . T x0 xk kb ON‐1 kt yk O0 Y0 k z0 zk Hình 4. Mô hình hệ vật – lò xo Sử dụng phương pháp Lagrange có thể đưa hệ phương trình chuyển động về dạng:     M(Φ, Φ, t )Φ  H(Φ, Φ, t )  Q(Φ, Φ, t ) (3) Phương pháp xoáy không dừng UVLM mở rộng được sử dụng để xác định lực khí động tác dụng lên cánh. Cánh được chia thành các panel có các khung xoáy tứ giác hoặc tam giác với lưu tốc bằng nhau. Tại các điểm kiểm tra trên mỗi panel, sử dụng điều kiện biên chảy không thấm Neumann để đưa ra hệ phương trình đại số nhằm xác định lưu tốc của mỗi khung xoáy. Phương pháp này đã được các tác giả sử dụng cho một số nghiên cứu trước đây cho cánh vẫy kiểu côn trùng [16-18] và cho kết quả tương đối chính xác. Hệ phương trình chuyển động (3) được giải bằng phương pháp Quasi-Newton dựa trên lý thuyết tìm đường do Tang và cộng sự đề xuất [19]. Mô hình giải thuật của chương trình tính toán được biểu diễn trên hình 5.
216 Nghiên cứu ảnh hưởng của độ cứng uốn đến lực nâng của cánh vẫy kiểu côn trùng Xây dựng mô hình kết cấu Mô hình FEM Xác định các đặc tính quán tính và độ cứng Tối ưu độ cứng Các thông số động học ban đầu Bộ giải Bộ giải khí động động lực học Tính trường áp suất Tải khí động trên bề mặt cánh Cập nhật Mô phỏng ĐLH đa vật mô hình khí động Điều kiện Không Trạng thái biến dạng hội tụ mới của cánh Thỏa mãn Trạng thái biến dạng của cánh, tải khí động và công suất t=t+dt Xuất dữ liệu Hình 5. Mô hình giải thuật 3. Kết quả Tiến hành khảo sát ảnh hưởng của độ cứng uốn đến lực nâng với các độ cứng xoắn khác nhau bằng cách thay đổi các hệ số 𝑎 và 𝑎 trong công thức (1). Ký hiệu 𝑎 và 𝑎 (𝑖 1, 2) là tỷ lệ giữa giá trị mới của 𝑎 và 𝑎 so với giá trị tương ứng của cánh ban đầu được cho trong bảng 1. Cho 𝑎 biến thiên từ 0,4 đến 2, 𝑎 thay đổi từ 0,9 đến 2,5. Cần chú ý, chỉ khảo sát 𝑎 từ 0,9 do ở các độ cứng xoắn nhỏ hơn, biến dạng xoắn của cánh trở nên quá lớn, không giống với kết cấu thực tế, cũng như lực pháp tuyến quá nhỏ nên không nằm trong phạm vi của bài toán thiết kế. Với mỗi giá trị của 𝑎 và 𝑎 xác định tần số dao động riêng thứ nhất f1 và lực nâng trung bình L của cánh trong một chu kỳ vẫy. Sự phụ thuộc của tần số dao động riêng thứ nhất f1 theo 𝑎 được biểu diễn trên hình 6. Tại mỗi giá trị của 𝑎 thể hiện cột sai số quanh giá trị trung bình, tương ứng với sự biến thiên của f1 theo các giá trị của 𝑎 . Khi độ cứng uốn tăng thì f1 tăng lên và sai số cũng tăng lên. Tuy nhiên có thể nói f1 chủ yếu phụ thuộc vào độ cứng uốn.
217 Lê Vũ Đan Thanh, Đặng Ngọc Thanh, Nguyễn Anh Tuấn 80 70 60 50 40 30 0.4 0.8 1.2 1.6 2 Hình 6. Tần số dao động riêng thứ nhất f1 của cánh với các độ cứng uốn khác nhau Trên hình 7 biểu diễn phân bố lực nâng trung bình trong một chu kỳ vẫy theo f1 với các giá trị khác nhau của độ cứng xoắn 𝑎 . Hình 7. Lực nâng trung bình theo f1 với các độ cứng xoắn khác nhau Có thể thấy, xu hướng biến thiên của lực nâng theo f1 tương đối giống nhau với các độ cứng xoắn khác nhau. Lực nâng tỉ lệ thuận với f1 khi độ cứng uốn nhỏ và ngược lại khi độ cứng uốn lớn. Sự phân bố này có thể được giải thích là do, khi độ cứng uốn giảm thì thành phần vận tốc theo phương ngang Vx tăng lên do biến dạng uốn. Khi đó vận tốc tại mút cánh tăng lên làm tăng lực pháp tuyến. Tuy nhiên, khi độ cứng uốn quá nhỏ sẽ khiến cánh có xu hướng bị uốn lên, tăng thành phần vận tốc theo phương thẳng đứng Vz. Khi đó, góc tấn hiệu dụng αe, được xác định theo công thức:  Vz   e   g  tan 1  , (4)  Vx  với αe là góc tấn hình học, sẽ bị giảm đi. Điều này đồng thời cũng làm giảm lực nâng. Các phân bố lực nâng theo f1 trong hình 7 được xấp xỉ bằng các đa thức bậc 3, và xác định giá trị 𝑓 tương ứng với độ cứng uốn tối ưu về lực nâng. Từ đó xác định tỉ lệ tần số tối ưu, là tỉ lệ giữa 𝑓 và tần số vẫy f, đối với mỗi giá trị độ cứng xoắn. Kết quả được biểu diễn trên hình 8. Có thể thấy, tỉ lệ tần số tối ưu với tất cả độ cứng xoắn đều xấp xỉ giá trị 2. Giá trị này gần với giá trị của loài nhậy Manduca Sexta trong tự nhiên, có tần số dao động riêng thứ nhất xấp xỉ 60 Hz và tần số vẫy 26,1 Hz.
218 Nghiên cứu ảnh hưởng của độ cứng uốn đến lực nâng của cánh vẫy kiểu côn trùng Nghiên cứu thực nghiệm của Ha [20] cho một số loài côn trùng tự nhiên như chuồn chuồn, ve sầu và bọ cánh cứng cũng cho thấy, lực nâng lớn nhất đạt được khi f1 xấp xỉ hai lần tần số vẫy. Hình 8. Tỉ lệ tần số tối ưu theo biên độ góc xoắn 4. Kết luận Trong bài báo này, các tác giả đã xây dựng sử dụng cách tiếp cận động lực học đa vật để xây dựng chương trình tính toán tương tác kết cấu – chất lưu cho cánh vẫy kiểu côn trùng. Cánh được mô hình hóa dưới dạng hệ vật – lò xo, kết hợp với tính toán khí động bằng phương pháp xoáy không dừng UVLM mở rộng. Chương trình được áp dụng để khảo sát ảnh hưởng của độ cứng uốn đến lực nâng trung bình của cánh. Khi độ cứng uốn của cánh nhỏ, lực nâng tỉ lệ thuận với độ cứng uốn và ngược lại ở độ cứng uốn lớn. Lực nâng đạt giá trị lớn nhất khi tần số dao động riêng thứ nhất xấp xỉ hai lần tần số vẫy. Kết quả này phù hợp với tỉ lệ tần số tối ưu của các loài côn trùng trong tự nhiên. Trong quá trình thiết kế TBB cánh vẫy phỏng côn trùng, cần tiến hành phân tích các dạng dao động riêng của kết cấu cánh, và lựa chọn tần số dao động riêng dựa trên tỉ lệ tần số tối ưu. Tài liệu tham khảo 1. Shyy, W., Aono, H., Kang, C.-k., and Liu, H. An Introduction to Flapping Wing Aerodynamics. Cambridge: Cambridge University Press, 2013. 2. Reid, H., Zhou, H., Maxcer, M., Peterson, R. K. D., Deng, J., and Jankauski, M. "Toward the design of dynamically similar artificial insect wings," International Journal of Micro Air Vehicles Vol. 13, 2021. doi: 10.1177/1756829321992138 3. Zhao, L., and Deng, X. "Power distribution in the hovering flight of the hawk moth Manduca sexta," Bioinspir Biomim Vol. 4, No. 4, 2009, p. 046003. doi: 10.1088/1748-3182/4/4/046003 4. Mountcastle, A. M., and Combes, S. A. "Wing flexibility enhances load-lifting capacity in bumblebees," Proc Biol Sci Vol. 280, No. 1759, 2013, p. 20130531. doi: 10.1098/rspb.2013.0531 5. Campos D, U. L., Bernal L. "Flow around flapping flexible flat plate wings," 50th Aiaa Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace
219 Lê Vũ Đan Thanh, Đặng Ngọc Thanh, Nguyễn Anh Tuấn Exposition. 2012. 6. Fu, J., Liu, X., Shyy, W., and Qiu, H. "Effects of flexibility and aspect ratio on the aerodynamic performance of flapping wings," Bioinspir Biomim Vol. 13, No. 3, 2018, p. 036001. doi: 10.1088/1748-3190/aaaac1 7. Du, G., and Sun, M. "Effects of wing deformation on aerodynamic forces in hovering hoverflies," J Exp Biol Vol. 213, No. Pt 13, 2010, pp. 2273-83. doi: 10.1242/jeb.040295 8. Truong, H., Engels, T., Kolomenskiy, D., and Schneider, K. "Influence of wing flexibility on the aerodynamic performance of a tethered flapping bumblebee," Theoretical and Applied Mechanics Letters Vol. 10, No. 6, 2020, pp. 382-389. doi: 10.1016/j.taml.2020.01.056 9. Le, V. D. T., Nguyen, A. T., Dang, N. T., and Phung, B. V. "A Multibody Dynamics Approach to Study an Insect-Wing Structure," Modern Mechanics and Applications. 2022, pp. 149-157. 10. Le, V. D. T., Nguyen, A. T., Phung, V. B., Pham, T. D., and Han, J.-H. "The beam modelling of the hawkmoth wing structure," Bioinspiration, Biomimetics, and Bioreplication X. 2020. 11. Le, V. D. T., Nguyen, A. T., and Dang, N. T. "Modal analysis of a body-spring system modeling an insect flapping-wing structure," Hội nghị khoa học các nhà nghiên cứu trẻ lần thứ XVI năm 2021. Học viện Kỹ thuật Quân sự, 2021. 12. O'Hara, R. P., and Palazotto, A. N. "The morphological characterization of the forewing of the Manduca sexta species for the application of biomimetic flapping wing micro air vehicles," Bioinspir Biomim Vol. 7, No. 4, 2012, p. 046011. doi: 10.1088/1748-3182/7/4/046011 13. Combes, S. A., and Daniel, T. L. "Flexural stiffness in insect wings. II. Spatial distribution and dynamic wing bending," J Exp Biol Vol. 206, No. Pt 17, 2003, pp. 2989-97. doi: 10.1242/jeb.00524 14. Pastor, M., Binda, M., and Harčarik, T. "Modal Assurance Criterion," Procedia Engineering Vol. 48, 2012, pp. 543-548. doi: 10.1016/j.proeng.2012.09.551 15. Willmott, A. P., and Ellington, C. P. "The mechanics of flight in the hawkmoth Manduca sexta. I. Kinematics of hovering and forward flight," Journal of Experimental Biology Vol. 200, No. 21, 1997, pp. 2705-2722. doi: 10.1242/jeb.200.21.2705 16. Nguyen, A. T., Kim, J.-K., Han, J.-S., and Han, J.-H. "Extended Unsteady Vortex- Lattice Method for Insect Flapping Wings," Journal of Aircraft Vol. 53, No. 6, 2016, pp. 1709-1718. doi: 10.2514/1.C033456 17. Nguyen, A. T., and Han, J.-H. "Wing flexibility effects on the flight performance of an insect-like flapping-wing micro-air vehicle," Aerospace Science and Technology Vol. 79, 2018, pp. 468-481. doi: 10.1016/j.ast.2018.06.007 18. Nguyen, A. T., Tran, N. D., Vu, T. T., Pham, T. D., Vu, Q. T., and Han, J.-H. "A Neural-network-based Approach to Study the Energy-optimal Hovering Wing
220 Nghiên cứu ảnh hưởng của độ cứng uốn đến lực nâng của cánh vẫy kiểu côn trùng Kinematics of a Bionic Hawkmoth Model," Journal of Bionic Engineering Vol. 16, No. 5, 2019, pp. 904-915. doi: 10.1007/s42235-019-0105-5 19. Tang, D., Bao, S., Lv, B., Guo, H., Luo, L., and Mao, J. "A derivative-free algorithm for nonlinear equations and its applications in multibody dynamics," Journal of Algorithms & Computational Technology Vol. 12, No. 1, 2017, pp. 30-42. doi: 10.1177/1748301817729990 20. Ha, N. S., Truong, Q. T., Goo, N. S., and Park, H. C. "Relationship between wingbeat frequency and resonant frequency of the wing in insects," Bioinspir Biomim Vol. 8, No. 4, 2013, p. 046008. doi: 10.1088/1748-3182/8/4/046008