Nghiên cứu đánh giá ảnh hưởng của một số thông số kết cấu cơ bản của khoang điều khiển đến ổn định bay của đạn phản lực

 lên từ 8o lên 11,25o tầm của đạn giảm đi cỡ 300m và nếu tăng từ 8o lên 14,5o tầm của đạn cỡ 700m. Kết quả này cho thấy góc mũi phần KĐK ảnh hưởng khá nhiều đến tầm bắn, do đó cần phải rất 0,36 14 ,5 o 11 ,2 5 o 8 o 0,2 0,26 120 chú ý trong quá trình thiết kế KĐK. Bảng 4.15 chỉ đưa ra các giá trị của hệ số tải khi  thay đổi, trong khi các thông số đánh giá khả năng điều khiển còn lại không được đưa ra do chúng thay đổi không đáng kể khi  thay đổi. Kết quả đưa ra ở bảng này cho thấy rằng khi  tăng lên vài độ thì hệ số tải giảm đi đáng kể, điều này có nghĩa là  ảnh hưởng tiêu cực đến khả năng điều khiển của đạn. Bảng 4.15. Hệ số tải nm thay đổi với nửa góc côn phần mũi KĐK  8o 11,25o 14,5o Hệ số tải nm 3,55 2,83 2,05 Hình 4.18. Thông số ổn định Ssd theo nửa góc mũi KĐK với góc lật cánh láiC = 10o Thông số ổn định động Ssd được thể hiện ở hình 4.18 khi  thay đổi. Các đường biểu diễn là khá gần nhau cho thấy sự thay đổi về tính chất ổn định của đạn là không nhiều, bảng 4.16 đưa ra thời điểm mất ổn định của các trường hợp tương ứng. Khi tăng  lên 6,5o (từ 8o lên 14,5o) thời điểm mất ổn định giảm đi 1s. Kết quả này cho thấy sự phụ thuộc yếu của góc côn phần mũi KĐK đến tính ổn định bay của đạn. Điều này cũng có thể giải thích rằng do thời điểm điều khiển được thực hiện khi đạn đã vượt qua đỉnh quỹ đạo và có vận tốc dưới âm. Tuy nhiên, nếu thay đổi thời điểm điều khiển sớm hơn khi đạn đang ở tốc độ trên âm thì đặc tính ổn định cũng như điều khiển của đạn sẽ 121 thay đổi rất nhiều vì ở vùng này tác động của góc  đến CD khá lớn như đã thể hiện ở bảng 4.13. Bảng 4.16. Thời điểm mất ổn định với nửa góc côn phần mũi KĐK  8o 11,25o 14,5o Thời điểm mất ổn định (s) 59,76 59,36 58,77 Nhận xét: Góc côn phần mũi KĐK càng nhỏ đạn càng ổn định, nhưng mức độ ảnh hưởng không nhiều, nhưng ảnh hưởng khá tiêu cực đến khả năng điều khiển của đạn. Kết luận chương 4 Qua các nội dung đã thực hiện ở ở chương 4, có thể rút ra một số kết luận sau: - Mô hình hình cơ sở của đạn phản lực 122mm đưa ra trong luận án có tính điều khiển phù hợp với loại đạn hiệu chỉnh quỹ đạo dựa trên cơ sở tính toán các thông số đánh giá khả khả năng điều khiển được đưa ra ở bảng 4.1. Do đó, có thể nghiên cứu ảnh hưởng một số thông số kết cấu cơ bản của KĐK đến ổn định bay của đạn phản lực 122mm lắp KĐK dựa trên mô hình này. - Biên độ góc lật cánh lái càng lớn đạn càng dễ mất ổn định và giá trị giới hạn của biên độ cánh được xác định qua công thức (4.7). Trong giới hạn của mô hình toán, cho thấy đạn mất ổn định khi điều khiển hướng còn điều khiển tầm thì rất ít ảnh hưởng. Khi điều khiển hướng hiệu chỉnh quỹ đạo với biên độ góc lật cánh tham chiếu C = 10o, không nên thực hiện quá 24,76s và nên thực hiện hiệu chỉnh quỹ đạo trong các khoảng thời gian ngắn với biên độ góc lật cánh lớn hơn tăng hiệu quả điều khiển mà đạn vẫn ổn định. - Cánh lái đặt càng gần xa tâm khối, đạn càng giảm ổn định và ngược lại. Khi cho vị trí cánh lái thay đổi 0 085% 105%C C Cx x x  (tương đương với thay đổi trong khoảng 0,2m) thì thời điểm mất ổn định thay đổi tăng 58,59s đến 63,38s, quy luật thay đổi khá tuyến tính khi cứ giảm 5% xC thì thời điểm mất 122 ổn định tăng cỡ 1s; độ dự trữ ổn định tĩnh SM tăng từ khoảng 8,36% ÷ 11,53% lên 9,9% ÷12,84% cho thấy đạn điều khiển tốn năng lượng hơn. - Diện tích cánh lái càng tăng, đạn càng dễ mất ổn định và ngược lại. Khi cho vị trí cánh lái thay đổi 20 085% 105%C C CS S S mm  thì thời điểm mất ổn định thay đổi giảm từ 63,68s đến 58,50s, quy luật thay đổi cũng khá tuyến tính khi cứ tăng 5% SC thì thời điểm mất ổn định tăng cỡ 1s; khoảng thay đổi độ dự trữ ổn định tĩnh SM là lớn khá lớn từ 1 1,8% (tương đương 10,28% ÷13,20% xuống 8,26% ÷ 11,42%) cho thấy tác động lớn đến năng lượng khi điều khiển. - Khối lượng của KĐK tăng lên sẽ làm đạn ổn định hơn, tuy khả năng điều khiển cũng như tầm của đạn bị giảm. Quy luật tăng của thời điểm mất ổn định so với khối lượng KĐK là khá phi tuyến, khi mF < 10 kg thì ảnh hưởng đến ổn định bay của đạn nhiều hơn khi mF > 10 kg, điều này cho thấy khối lượng KĐK mF = 10kg là giá trị hợp lý. - Tâm khối KĐK càng gần mũi đạn đạn càng ổn định hơn, nhưng mức ảnh hưởng là không thực sự rõ rệt, khi khảo sát vị trí tâm khối KĐK từ 0,256m đến 0,456m thì thời điểm mất ổn định giảm từ 60,38s đến 59,13s, tuy nhiên lại ảnh hưởng khá nhiều đến độ ổn tĩnh với lượng thay đổi cỡ 1,5% khi thay đổi tâm khối KĐK trong giới hạn 0,2m. Điều này cho thấy rằng khi thiết kế cần phải bố trí khối tâm của KĐK gần mũi nhất có thể để giảm năng lượng điều khiển cánh lái. - Góc côn phần mũi KĐK không ảnh hưởng nhiều đến ổn định của đạn nhưng lại ảnh hưởng nhiều đến khả năng điều khiển đạn ở quỹ đạo có bán kính nhỏ thể hiện khi tăng nửa góc côn phần mũi đạn từ 8o lên 14,5o thời điểm mất ổn định giảm từ 59,76s xuống 58,77s và hệ số tải nm giảm từ 3,55 xuống 2,05. Đối với đạn hiệu chỉnh quỹ đạo lượng hiệu chỉnh thường là nhỏ do đó việc tăng góc côn phần mũi vẫn được cho phép nếu không gian bố trí bên trong cần thiết. 123 Nghiên cứu ảnh hưởng của một số thông số cơ bản của KĐK đến ổn định bay và khả năng điều khiển của đạn phản lực 122mm lắp KĐK, đã đưa ra mức độ ảnh hưởng của các thông số trong giới hạn cho phép của các thông số phù hợp với đối tượng, từ đó cho thấy ý nghĩa của kết quả nghiên cứu cho quá trình thiết kế KĐK nói riêng và đến cả quả đạn nói chung. 124 KẾT LUẬN Kết quả nghiên cứu của luận án: Từ các kết quả nghiên cứu, luận án rút ra một số kết luận sau: 1. Đạn phản lực được lắp KĐK với đặc điểm là cơ hệ 7 bậc tự do có đặc tính ổn định phức tạp, do đó nghiên cứu ổn định bay của loại đạn này là rất cần thiết để phục vụ cải tiến tăng độ chính xác cho loại đạn này. 2. Mô hình toán mô tả chuyển động bay của đạn phản lực lắp KĐK được xây dựng trong luận án bao gồm 14 phương trình vi phân độc lập và 14 biến trạng thái độc lập được thể hiện ở (2.35), (2.38), (2.46) và (2.47) được xây dựng trên cơ sở khoa học chắc chắn và mô tả được chuyển động phức tạp 7 bậc tự do của đạn trên quỹ đạo. Mô hình này có thể dùng cho nhiều đối tượng đạn lắp KĐK khác nhau. 3. Điều kiện ổn định bay của đạn phản lực lắp KĐK được xây dựng trong luận án bao gồm 2 bất đẳng thức (2.98) và (2.100) là phù hợp với đối tượng đạn 7 bậc tự do của luận án. Điều kiện ổn định bay này có tính tổng quát cao vì được thiết lập ở góc bắn bất kỳ và không bỏ qua tác động của trọng lực như đối với điều kiện ổn định bay của đạn 6 bậc tự do truyền thống, do đó điều kiện ổn định này có thể áp dụng cho nhiều đối tượng đạn khác nhau. 4. Luận án đã chứng minh được tính tin cậy của mô hình toán mô tả chuyển động bay và điều kiện ổn định bay của đạn phản lực lắp KĐK bằng các luận giải có căn cứ khoa học và so sánh với các tài liệu được công bố trên các tạp chí khoa học uy tín. 5. Nghiên cứu được thực hiện ở chương 4 của luận án đã cho thấy ảnh hưởng của một số thông số kết cấu cơ bản của KĐK đến ổn định bay của đạn là khác nhau ở mức độ, chiều hướng tác động cụ thể là: - Biên độ góc lật cánh lái, khoảng cách cánh lái đến tâm khối và diện tích cánh lái càng lớn đạn càng dễ mất ổn định và ngược lại. Diện tích cánh lái ảnh 125 hưởng nhiều tới ổn định bay và khả năng điều khiển của đạn khi tác động tới cả hệ số lực và mô men khí động cánh lái. - Khối lượng KĐK tăng hoặc tâm khối KĐK càng gần mũi đạn đạn càng ổn định và ngược lại. Tâm khối KĐK ảnh hưởng không nhiều đến ổn định nhưng lại ảnh hưởng khá nhiều đến độ dự trữ ổn định tĩnh, do đó nên bố trí tâm khối KĐK gần mũi đạn nhất có thể. - Góc côn phần mũi KĐK không ảnh hưởng nhiều đến ổn định của đạn nhưng lại ảnh hưởng nhiều đến khả năng điều khiển đạn ở quỹ đạo có bán kính nhỏ. Đối với đạn hiệu chỉnh quỹ đạo lượng hiệu chỉnh thường là nhỏ do đó việc tăng góc côn phần mũi vẫn được phép nếu không gian bố trí bên trong cần thiết. Những đóng góp mới của luận án: 1. Xây dựng mô hình toán mô tả chuyển động 7 bậc tự do trên quỹ đạo của đạn phản lực lắp khoang điều khiển, trên cơ sở đó xây dựng điều kiện ổn định bay cho mô hình nghiên cứu. 2. Đánh giá ảnh hưởng của một số thông số kết cấu cơ bản của khoang điều khiển đến ổn định bay của đạn phản lực lắp khoang điều khiển. Hướng nghiên cứu tiếp theo: 1. Mở rộng mô hình toán và xây dựng điều kiện ổn định bay cho đạn phản lực ổn định quay lắp KĐK. 2. Nghiên cứu xây dựng bộ điều khiển hiệu chỉnh quỹ đạo dựa trên mô hình toán đã xây dựng. Nghiên cứu xác định thời điểm điều khiển hợp lý cho đạn trên quỹ đạo. 3. Tối ưu hóa thông số kết cấu KĐK và nghiên cứu bổ sung một số thông số kết cấu của KĐK như hình dạng cánh lái, chiều dày cánh lái, chiều dài KĐK... tác động đến ổn định bay của đạn. 126 DANH MỤC CÔNG TRÌNH KHOA HỌC ĐÃ CÔNG BỐ 1. Trần Xuân Diệu, Nguyễn Phú Thắng, Phan Văn Chương, Vũ Minh Đức, "Mô hình hóa và mô phỏng quỹ đạo bay của đạn pháo phản lực có ngòi hiệu chỉnh quỹ đạo dạng tách chuyển động quay," TC. Nghiên cứu KHCNQS, số 54, tr. 22-32, Hà Nội, 4-2018. 2. Trần Xuân Diệu, Nguyễn Phú Thắng, Phan Văn Chương, Trần Quang Minh, "Xây dựng điều kiện ổn định bay của đạn pháo phản lực," TC. Nghiên cứu KHCNQS, số đặc san, tr. 205-214, Hà Nội, 8-2018. 3. Trần Xuân Diệu, Nguyễn Phú Thắng, Phan Văn Chương, "Xây dựng điều kiện ổn định bay của đạn pháo phản lực lắp khoang điều khiển," TC. Nghiên cứu KHCNQS, số đặc san, tr. 232-241, Hà Nội, 10-2019. 4. Trần Xuân Diệu, Nguyễn Phú Thắng, Phan Văn Chương, "Nghiên cứu ảnh hưởng của một số thông số kết cấu khoang điều khiển đến tính ổn định trên quỹ đạo của vật bay lắp khoang điều khiển bằng mô hình góc tấn phức," TC. Cơ khí Việt Nam, số 6, tr. 88-92, Hà Nội, 6-2020. 5. Trần Xuân Diệu, Nguyễn Phú Thắng, Phan Văn Chương, Nguyễn Mạnh Tuấn, "Nghiên cứu ảnh hưởng của góc lật cánh lái đến ổn định bay của đạn pháo phản lực lắp khoang điều khiển," TC. Nghiên cứu KHCNQS, số đặc san, tr. 262-267, 10-2020. 127 TÀI LIỆU THAM KHẢO Tiếng Việt 1. Bảng bắn pháo phản lực BM-21 (2003), Nhà xuất bản quân đội nhân dân, Hà Nội. 2. Bộ môn đạn, Khoa Trang bị cơ điện (1976), Nguyên lý kết cấu và tính toán đạn phản lực, Trường Đại học Kỹ thuật Quân sự. 3. Bộ môn Thuật phóng và Điều khiển hỏa lực, Khoa Vũ khí (2004), Cơ sở thiết kế tên lửa, Học viện Kỹ thuật Quân sự. 4. Bộ môn đạn, Khoa Vũ khí (2006), Cơ sở tính toán đạn phản lực không điều khiển, Học viện Kỹ thuật Quân sự. 5. Cục Quân khí (2005), Quy trình công nghệ kiểm nghiệm pháo phản lực 122mm. 6. Nguyễn Thanh Hải (2014), Nghiên cứu nâng cao độ ổn định của tên lửa không điều khiển bằng lựa chọn tham số cánh hợp lý, Luận án tiến sỹ kỹ thuật, Hà Nội. 7. Nguyễn Lạc Hồng, Nguyễn Xuân Anh (1998), Giáo trình tính toán bệ phóng và động lực học khi phóng, Học viện Kỹ thuật Quân sự, Hà Nội. 8. Nguyễn Xuân Hùng (2004), Động lực học và ổn định bay, Nhà XBQG, Hà Nội. 9. Nguyễn Văn Khang (2017), Động lực học hệ nhiều vật, Nhà xuất bản KHKT, Hà Nội. 10. Nguyễn Văn Khối (2016), Tìm hiểu thông tin về khả năng thông minh hóa cho đạn phản lực tăng tầm BM-21 và đạn pháo nòng dài tăng tầm 130mm M46, Viện Khoa học và Công nghệ Quân sự, Hà Nội. 11. Phạm Văn Luận (2013), Xác định thời điểm làm việc hợp lý của động cơ tên lửa nhiên liệu rắn cỡ nhỏ được phóng từ giàn phóng theo nguyên lý khí động cao áp, thấp áp, Luận án Tiến sỹ Kỹ thuật, Hà Nội. 12. Nguyễn Duy Phồn (2016), Nghiên cứu ảnh hưởng của một số yếu tố 128 nhiễu động ban đầu khi đạn tên lửa không điều khiển rời bệ phóng đến độ chính xác bắn, Luận án tiến sĩ kỹ thuật, Hà Nội. 13. Lê Quang (2016), Cơ học vật bay, Đại học Bách khoa Hà Nội, Hà Nội. 14. Trần Ngọc Quý (2018), "Phương pháp hiệu chỉnh quỹ đạo đạn pháo phản lực kiểu 9M22Y cải tiến", Tạp chí nghiên cứu khoa học và công nghệ quân sự, số 58, tr. 3-13. 15. Trần Ngọc Quý (2018), Nghiên cứu khả năng điều khiển của đạn pháo phản lực 9M22Y trên cơ sở ứng dụng hệ thống xử lý thông tin kết hợp INS/GPS(GLONASS), Đề tài nền cấp Viện KH-CN Quân sự, Hà Nội. 16. Lê Minh Thái (1999), Giải bài toán chuyển động của đạn khi kể đến các yếu tố nhiễu loạn, Luận án Tiến sỹ kỹ thuật, Hà Nội. 17. Nguyễn Văn Thọ (2003), Giáo trình thuật phóng ngoài, Học viện Kỹ thuật Quân sự, Hà Nội. 18. Nguyễn Đức Thuận (2018), Nghiên cứu, xây dựng phương án giảm tản mát cho đạn phản lực tăng tầm 40km bắn trên dàn phóng xe chiến đấu BM-21, Đề tài cấp Viện KH&CN QS, Hà Nội. 19. Viện Thiết kế Vũ khí (1995), Bản vẽ sản phẩm đạn phản lực 122mm, Tổng cục Công nghiệp Quốc phòng và Kinh tế. Tiếng Anh 20. J. D. Anderson (1999), Aircraft Performance and Design, McGraw-Hill, London. 21. J. S. Barrowman (1967), The practical calculation of the aerodynamic characteristics of slender finned vehicles, Catholic University of America, Washington. 22. M. E. Bartowitz (2008), Determination of Static and Dynamic Stability Coefficients Using Beggar, Air Force Institute of Technology, USA. 23. W. B. Blake (1998), Missile Darcom user's manual - 1997 fortran 90 revision, Air Force Research Laboratory, USA. 129 24. S. S. Chin (1961), Missile Configuration Design, McGraw-Hill, Orlando, Floria, USA. 25. M. V. Cook (2007), Flight Dynamics Principles, Butterworth- Heinemann, Burlington, USA. 26. M. Costello (2000), Modelling and simulation of a differential roll projectile, Army Research Laboratory, USA. 27. M. Costello, A. Peterson (2000), "Linear theory of a dual-spin projectile in atmospheric flight", Journal of guidance and control, vol. 23, no. 5, pp. 789-797. 28. O. Dullum (2010), The Rocket Artillery Reference Book, Norwegian Defence Research Establishment, USA. 29. B. Etkin, L. D. Reid (1995), Dynamics of Flight: Stability and Control, John Wiley and Sons, New York, USA. 30. J. Ferris (1967), Static stability investigation of a single- stage sounding rocket at Mach numbers from 0.60 to 1.20, NASA-TN-D-4013, USA. 31. E. L. Fleeman (2001), Tactical Missile Design, American Institute, Virginia, USA. 32. A. E. Gamble, P. N. Jenkins (2002), "Low cost guidance for the Multiple Launch Rocket System (MLRS) artillery rocket", IEEE Aerospace and Electronic Systems Magazine, vol. 6, no. 1, pp. 33-39. 33. M. Gopal (2002), Control Systems: Principles and Design, Tata McGraw-Hill Education, New Delhi. 34. Q. Guo (2015), "Guidance law design for a class of dual-spin mortars", International Journal of Aerospace Engineering, vol. 2015, pp. 1-12. 35. M. J. Hemsch (1992), Component Build-Up Method for Engineering Analysis of Missiles at Low-to-High Angle of Attack, Hampton, Virginia, USA. 36. D. E. Hoak et al (1960), USAF Stability an Control Datcom, Air Force 130 Wright Aeronautical Laboratories, USA. 37. L. H. Johson (2012), Precision Guidance Kit, 56th Annual Fuze Conference, Baltimore. 38. D. Karnopp (2004), Vehicle Stability, University of California, California, USA. 39. K. H. Lloyd, D. P. Brown (1977), Influence of gravity and applied side forces on the stability of a spinning projectile, Defence Science and Technology Oganisation, Australia. 40. K. H. Lloyd (1981), "Instability of spinning projectiles during terminal guidance", Journal of Guidance, Control, and Dynamics, vol. 1, no. 4, pp. 66-69. 41. M. Martinez (2016), Preliminary aeromechanical design of sub-guided munition for artillery rockets, 29th International Symposium on Ballistics, Edinburgh, Scotland. 42. R. L. McCoy (1999), Modern Exterior Ballistics, Schiffer Edition, Pennsylvania, USA. 43. Military Handbook (1990), Design of aerodynamically stabilized free rockets, Rocket Science Institute, USA. 44. H. Molitz, R. Strobel (1963), Äussere Ballistik, Springer Verlag, Berlin, German. 45. C. H. Murphy (1963), Free flight motion of symetric missiles, BRL, USA. 46. C. H. Murphy (1979), Effect of horizontal and vertical side forces and moments on stability of a symmetric missile in ascending or descending flight, Ballistic Research Laboratory, USA. 47. C. H. Murphy (1981), "Instability of controlled projectiles in ascending or descending flight", Journal of Guidance and Control, vol. 4, no. 1, pp. 66-69. 131 48. P. Myrick (2010), Affordable Precision Artillery Precision Guidance Kit (PGK), 8th Annual Future Artillery Conference, London. 49. J. D. Nicholaides (1959), Missile flight and astrodynamics, Bureau of Naval Weapons, Washington, USA. 50. W. Pettis (1967), Aerodynamic Characteristics of Multiple Fins of Rectangular Planform on a Body of Revolution at Mach Numbers of 1.48 to 2.22, US Army Missile Command, USA. 51. C. N. Robert (1998), Flight Stability and Automatic Control, McGraw- Hill, USA. 52. F. Shahid (2017), "Variation in aerodynamic coefficients with altitude", Results in Physics, vol. 7, pp. 1261-1273. 53. Standard atmosphere (1975), International organization for standardization, ISO 2533:1075. 54. S. Theodoulis et al. (2013), "Guidanceand control design for a class of spin-stabilized fin-controlled projectiles", Journal of Guidance, Control, and Dynamics, vol. 36, no. 2, pp. 517-531. 55. S. Theodoulis, P. Wernert (2011), Modelling and stability analysis for a class of 155mm spin-stabilized projectiles with course correction fuse, AIAA atmospheric flight mechanics conference and exhibit, Portland, Oregon, USA. 56. Y. Wang et al. (2015), "Correction mechanism analysis for a class of spin-stabilized projectile with fixed canards", Engineering Letters, vol. 23, no. 4, pp. 269-276. 57. Y. Wang et al. (2018), "Influence of control strategy on stability of dual- spin projectiles with fixed canards", Defence Technology, vol. 14, no. 6, pp. 709-719. 58. P. Wernert (2009), Stability analysis for canard guided dual-spin stabilized projectiles, AIAA atmospheric flight mechanics conference 132 and exhibit, Chicago, USA. 59. P. Wernert (2011), Modeling and stability analysis for a class of 155 mm spin-stabilized projectiles with course correction fuse(CCF), AIAA Atmospheric Flight Mechanics Conference and Exhibit, Oregon, USA. 60. G. L. Winchenbach (1997), Aerodynamic Testing in a Free-Flight Spark Range, Wright Laboratory Armament Directorate, USA. 61. D. Zhu et al. (2014), "Flight stabitity of a dual-spin projectile with canards", Journal of Aerospace Engineering, vol. 229, no. 4, pp. 703- 716. Tiếng Nga 62. Боевая машина БМ21, Технические условия на капитальный ремонт БМ-21 УК (1980), Военное издательство, Москва. 63. А. А. Дмитриевский (1972), Внешняя баллистика, Машиностроение, Москва. 64. О. В. Коротков (2016), Aвтономный блок управления реактивного снаряда, Федеральная служба по интеллектуальной собственности, Правительство России. 65. О. В. Коротков, С. Г. Благов (2016), "Дальнобойный высокоточный управляемый артиллерийский снаряд большого калибра с аппаратурой спутниковой радионавигации", Известия Российской Академии Ракетных и Артиллерийских Наук, № 1, Ч. 1, C. 60-66. 66. О. В. Коротков (2017), "Aвтономные отсеки управления с аппаратурой спутниковой навигации для 122-мм реактивных снарядов. Kомплексная математическая модель", Известия ТулГУ. Технические Hауки, № 4, Ч. 3, C. 135-145. 67. Л.М. Левин, Ф. Р. Гантмахер (1959), Теория полета неуправляемых ракет, Физматгиз, Москва. 68. С. С. Максимов, Д. Н. Автух (2012), "Aнализ возможных схем 133 построения систем управления рс повышенной точности для рсзо калибра 122-мм", Известия ТулГУ. Технические науки, № 11, Ч. 1, C. 167-170. 69. Я.М. Шапиро (1951), Пороховые реактивные снаряды, Артиллерийской Академии, Москва. 70. К. Яцек (1988), Летающие крылья, ДОСААФ СССР, Москва. P1 PHỤ LỤC 1 Chương trình tính toán thông số chuyển động của đạn trong ống phóng 1.1. Chương trình chính clc; clear all; close all; v0 = 0;%van toc ban dau xd0 = 0; %vi tri ban dau deltarx = 1.5*pi/180; d = 0.122; Ld = 3; %chieu dai OP x0 = [xd0 v0 ]; deltat = 0.0001; Tf = 0.5; [Td,Xd,Ts_xuat] = rungekutta(Tf,x0,deltat,@hptdan_chuyen_dong_trong_nong); x = Xd(:,1);%quang duong v = Xd(:,2);%van toc dan nd= find(x-Ld<0,1,'last'); %kich thuoc vec-tor khi dan con trong OP Td = Td(1:nd); x=x(1:nd); v=v(1:nd); omega_opA = 2*v/d*tan(deltarx); v0 = v(nd) omega0 = omega_opA(nd) figure(1); subplot(2,2,1); plot(Td,x);title('Quang duong'); grid on; subplot(2,2,2); plot(Td,v);title('Van toc dai');grid on; subplot(2,2,3); plot(Td,omega_opA);title('Toc do quay quanh truc dan');grid on; figure(2) plot(Td,v);grid on; xlabel('Thoi gian,[s]'); ylabel('Van toc cua dan trong OP [m/s]') figure(3) plot(Td,omega_opA); grid on; xlabel('Thoi gian,[s]'); ylabel('Toc do quay cua dan trong OP [rad/s]') 1.2. Hệ phương trình vi phân mô tả chuyển động đạn trong ống phóng function [dx,g] = hptdan_chuyen_dong_trong_nong(t,x) dx=zeros(2,1); tc = 1.88; %thoi gian chay cua thuoc phong mtp = 20.5; % khoi luong cua thuoc phong mf = 0;%10;%khoi luong phan truoc ma = 66.6;%khoi luong ban dau phan sau Fdc = 23000; %luc day dong co md = -mtp/tc*t + ma + mf; deltarx = 1.5*pi/180; Ixa = -0.0144*t+0.15; d = 0.122; g = 9.8; muyrx = 0.002; temp1 = cot(deltarx) - muyrx - (tan(deltarx)+muyrx)*4*Ixa/md/d^2; temp2 = cot(deltarx) - muyrx + (tan(deltarx)+muyrx)*4*Ixa/md/d^2; dx(1) = x(2); dx(2)=(temp1/temp2)*(Fdc/md); end 1.3. Thuật toán Rung-Kutta function [T,X,Ts_xuat] = rungekutta(Tf,X0,deltat,hpt) P2 T = 0:deltat:Tf; n = length(T); t = 0; x = X0'; m = length(X0); X = zeros(n,m); Ts_xuat = zeros(n,1); for i = 1:n for j = 1:m X(i,j) = x(j); end [k1,temp] = hpt(t,x); Ts_xuat(i) = temp; k2 = hpt(t+deltat/2,x+deltat/2*k1); k3 = hpt(t+deltat/2,x+deltat/2*k2); k4 = hpt(t+deltat,x+deltat*k3); x = x + deltat/6*(k1+2*k2+2*k3+k4); t = t +deltat; end PHỤ LỤC 2 Chương trình tính toán thông số chuyển động bay của đạn phản 122mm 2.1. Chương trình chính clc; clear all;% close all; u0 = 45.82;%Duoc giai tu phuong trinh vi phan trong OP v0 = 0; w0 = 0; theta0 = 45.78*pi/180; % goc ban ban dau psi0 = 0; %goc huong ban dau xe0 = 0; ye0 = 0; ze0 = 0; p0 = 19.67; %rad/s toc do quay ban dau. Duoc giai tu phuong trinh vi phan trong OP q0 = 0; % toc do goc pitch ban dau r0 = 0; %toc do goc yaw ban dau x0 = [u0 v0 w0 theta0 psi0 xe0 ye0 ze0 p0 q0 r0]; %index 1 2 3 4 5 6 7 8 9 10 11 global deltat deltat = 0.01; Tf = 85; [T,X,Ts_xuat] = rungekutta(Tf,x0,deltat,@bm21nr); u = X(:,1); v = X(:,2); w = X(:,3); theta = X(:,4); psi = X(:,5); xe = X(:,6); ye = X(:,7); ze = X(:,8); p = X(:,9); q = X(:,10); r = X(:,11); V = sqrt(u.^2+v.^2+w.^2); alpha = w./V; beta = v./V; alphat = sqrt(alpha.^2+beta.^2); nd= find(-ze>0,1,'last'); %kich thuoc vec-tor khi dan cham muc tieu tam_ban = xe(nd) % tam ban lon nhat thoi_gian_bay = T(nd) P3 van_toc_roi = V(nd) goc_roi = theta(nd)*180/pi do_cao_duong_dan = max(-ze(1:nd)) van_toc_lon_nhat = max(V) T = T(1:nd); % tong thoi gian dan bay %----------------------------------------------------------------------- xz= [xe(1:nd),-ze(1:nd)]; savefile = 'xz4578.mat'; save(savefile,'xz'); %------------------------------------------------------------------------ subplot(2,2,1); plot(T,alpha(1:nd)*180/pi,'Color','black','Linewidth',1.5);grid on; ylabel('goc tan, alpha[do]');xlabel('thoi gian, t[s]'); hold on; subplot(2,2,2); plot(T,beta(1:nd)*180/pi,'Color','black','Linewidth',1.5);grid on; ylabel('goc truot canh, beta[do]');xlabel('thoi gian, t[s]'); hold on; subplot(2,2,3); plot(xe(1:nd),-ze(1:nd));xlabel('Tam');ylabel('Cao do');title('Mat phang ban');grid on; subplot(2,2,4); plot3(xe(1:nd),ye(1:nd),-ze(1:nd)); xlabel('Tam');ylabel('Huong');zlabel('Cao do');title('Khong gian 3 chieu'); grid on 2.2. Hệ phương trình vi phân chuyển động bay đạn phản lực 122mm function [dx,g] = bm21nr(t,x) dx = zeros(11,1); global deltat m = (-11.4*t +66.6)*(t=1.88); g = 9.8; Ix = (-0.0144*t+0.15)*(t=1.74); Iy = (-4.3*t+41.58)*(t=1.74); d = 0.122; l = d; S = 1/4*pi*d^2; Maro = khiquyen(-x(8)); % Toc do am thanh va mat do theo do cao ro = Maro(1); at = Maro(2); %vt am thanh V = sqrt(x(1)^2+x(2)^2+x(3)^2); M = V/at; %hs Mach alpha = x(3)/x(1); beta = x(2)/V; %----------Hs khi dong ----------------- Adnm = hskd9M22U(M); CD = 1*Adnm(1); CLalpha = Adnm(2); Cl = Adnm(3); Clp = Adnm(4); CMalpha = Adnm(5); CMq = Adnm(6); CMalpha_dot = Adnm(7); hsl = 1/2*ro*S*V^2; hsm = hsl*l; L = hsl*CLalpha*[0;-beta;-alpha;];% luc nang tong hop D = -hsl*CD*[1; 0; 0;]; % luc can tong hop[1; beta; alpha;] K = [0;0;0;];%hsl*Cypalpha*(x(9)*l/V)*[0;-alpha; beta;];% luc Magnus tong hop Fdc = luc_day_dc(t)*[1;0;0;]; Mp = hsm*CMalpha*[0;alpha;-beta;]; % momen lat tong hop Mm = [0;0;0;];%hsm*(-CMpalpha*x(9)*l/V)*[0;-beta;-alpha;];% mo men magnus tong hop Mpd = hsm*(CMq + CMalpha_dot)*l/V*[0;x(10);x(11);]; % momen giam chan tong hop Mrd = hsm*Clp*x(9)*l/V*[1;0;0;]; % MM giam chan roll cho phan sau Ml = hsm*Cl*[1;0;0;]; %MM can roll tong quat P4 Ft = L + D + K + Fdc; Mt = Mp + Mm + Mpd + Mrd + Ml; X = Ft(1); Y = Ft(2); Z = Ft(3); Mx = Mt(1); My = Mt(2); Mz = Mt(3); dx(1) = X / m - g * sin(x(4)) + x(11) * x(2) - x(10) * x(3); dx(2) = Y / m - x(11) * x(1) - x(11) * tan(x(4)) * x(3); dx(3) = Z / m + g * cos(x(4)) + x(10) * x(1) + x(11) * tan(x(4)) * x(2); dx(4) = x(10); dx(5) = 0.1e1 / cos(x(4)) * x(11); dx(6) = cos(x(4)) * cos(x(5)) * x(1) - sin(x(5)) * x(2) + sin(x(4)) * cos(x(5)) * x(3); dx(7) = cos(x(4)) * sin(x(5)) * x(1) + cos(x(5)) * x(2) + sin(x(4)) * sin(x(5)) * x(3); dx(8) = -sin(x(4)) * x(1) + cos(x(4)) * x(3); dx(9) = 0.1e1 / Ix * Mx; dx(10) = 0.1e1 / Iy * (My - x(11) * Ix * x(9) - x(11) ^ 2 * tan(x(4)) * Iy); dx(11) = 0.1e1 / Iy * (Mz + x(10) * Ix * x(9) + x(11) * tan(x(4)) * Iy * x(10)); 2.3. Hệ số khí động đạn phản lực 122mm function Adnm = hskd9M22U(M) Adnm=zeros(1,7); a = [0.2 0.378 8.57 0.075 -4.099 -53.8386 -1271.46 -48.2965 0.4 0.366 9.1 0.076 -4.156 -54.3829 -1273.345 -49.7215 0.6 0.362 9.48 0.077 -4.2245 -55.1018 -1269.6 -52.2935 0.8 0.376 9.82 0.076 -4.224 -54.3483 -1236.84 -56.314 1 0.505 9.99 0.077 -4.1695 -55.0595 -1211.535 -61.999 1.1 0.579 10.16 0.094 -4.337 -62.072 -1337.76 -65.6125 1.2 0.578 10 0.102 -4.435 -65.6489 -1390.365 -70.0445 1.3 0.562 10.08 0.11 -4.652 -66.055 -1438.53 -72.24 1.4 0.546 10.34 0.123 -5.085 -73.4021 -1534.215 -73.7205 1.5 0.530 10.44 0.128 -4.9735 -70.0033 -1140.315 -74.4035 1.6 0.514 10.59 0.126 -4.918 -68.1121 -1141.22 -74.8755 1.8 0.484 10.67 0.12 -4.6695 -63.0248 -1087.59 -75.431 2 0.457 10.13 0.113 -4.381 -57.607 -1094.02 -75.6045 2.2 0.432 9.05 0.105 -4.075 -37.6782 -920.41 - 75.4765;]; if M<a(1,1) for j = 2:8 m = a(1,j)/a(1,1); Adnm(j-1)=m*M; end end if M>a(14,1) Adnm = a(14,2:8); end if M<=a(14,1) for i = 1:14 if M==a(i,1) Adnm = a(i,2:8); elseif M>a(i,1)&&M<a(i+1,1) for j = 2:8 k = (a(i+1,j)-a(i,j))/(a(i+1,1)-a(i,1)); Adnm(j-1)=k*(M-a(i,1))+a(i,j); end end end P5 end 2.4. Mật độ và tốc độ âm thanh theo mô hình khí tượng tiêu chuẩn 1975 function Maro = khiquyen(M) Maro=zeros(1,2); % Do cao Mat do Toc do am thanh a = [ 0 1.225 340.29 500 1.1673 338.37 1000 1.1117 336.44 1500 1.0581 334.49 2000 1.0066 332.53 2500 0.95695 330.56 3000 0.90925 328.58 3500 0.8634 326.59 4000 0.81935 324.59 4500 0.77704 322.57 5000 0.73643 320.55 6000 0.66011 316.45 7000 0.59002 312.31 8000 0.52579 308.11 9000 0.46706 303.85 10000 0.41351 299.53;]; if M<a(1,1) for j = 2:3 m = a(1,j)/a(1,1); Maro(j-1)=m*M; end end if M>a(16,1) Maro = a(16,2:3); end if M<=a(16,1) for i = 1:16 if M==a(i,1) Maro = a(i,2:3); elseif M>a(i,1)&&M<a(i+1,1) for j = 2:3 k = (a(i+1,j)-a(i,j))/(a(i+1,1)-a(i,1)); Maro(j-1)=k*(M-a(i,1))+a(i,j); end end end end PHỤ LỤC 3 Chương trình tính toán thông số chuyển động bay của đạn pháo 155mm 3.1. Chương trình chính clc; clear all; close all; u0 = 500; v0 = 0; w0 = 0; phif0 = 0; % goc quay ban dau cua phan phia truoc phia0 = 0; % goc quay ban dau cua phan phia sau theta0 = 40*pi/180; % goc ban ban dau psi0 = 0; %goc huong ban dau xe0 = 0; ye0 = 0; ze0 = 0; P6 pf0 = 1445; % rad/s toc do ban dau phan trc pa0 = 1445; %rad/s toc do ban dau phan sau q0 = 0; % toc do goc pitch ban dau r0 = -1.5; %toc do goc yaw ban dau x0 = [u0 v0 w0 phif0 phia0 theta0 psi0 xe0 ye0 ze0 pf0 pa0 q0 r0 0 0 0 0 0 0]; %index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 deltat = 0.01; Tf = 50; [T,X,Ts_xuat] = rungekutta(Tf,x0,deltat,@hptdualspin155); u = X(:,1); v = X(:,2); w = X(:,3); phif = X(:,4); phia = X(:,5); theta = X(:,6); psi = X(:,7); xe = X(:,8); ye = X(:,9); ze = X(:,10); pf = X(:,11); pa = X(:,12); q = X(:,13); r = X(:,14); FCx = X(:,15); FCy = X(:,16); FCz = X(:,17); MCx = X(:,18); MCy = X(:,19); MCz = X(:,20); V = sqrt(u.^2+v.^2+w.^2); alpha = w./V; beta = v./V; alphat = sqrt(alpha.^2+beta.^2); subplot(2,2,1); plot(T,alpha*180/pi,T,alphat*180/pi);title('alpha-goc tan'); subplot(2,2,2); plot(T,beta*180/pi);title('beta-goc truot canh'); subplot(2,2,3); plot(xe,-ze);title('ye'); subplot(2,2,4); plot3(xe,ye,-ze); xlabel('Tam');ylabel('Huong');zlabel('Cao do') grid on 3.2. Hệ phương trình vi phân chuyển động bay đạn pháo 155mm function [dx,g] = hptdualspin155(t,x) dx = zeros(17,1); mf = 1.482; ma = 41.316; m = mf + ma; Ixf = 0.014; Ixa = 0.133; Iy = 1.893; l = 0.155; S = 0.01887; g = 9.8; Cv = 0.001; CR = 0.001; ro = 1.2; at = 330; %vt am thanh V = sqrt(x(1)^2+x(2)^2+x(3)^2); M = V/at; %hs Mach alpha = x(3)/x(1); beta = x(2)/V; %----------Hs khi dong ----------------- P7 CM155 = hskddual155(M); CLalpha = CM155(3); CD = CM155(1)+CM155(2); CMalpha = CM155(4); CMq = CM155(5); CMalpha_dot = CM155(6); Clp_F = CM155(7); Clp_A = CM155(8); Cypalpha_F = CM155(9); Cypalpha_A = CM155(10); CMpalpha_F = CM155(11); CMpalpha_A = CM155(12); CNdelta = CM155(13); CMdelta = CM155(14); %---------Goc lat canh lai ----------------- phiC = pi/2; delta = 12*pi/180*(t>25); % deltay = delta*cos(rem(x(4),2*pi)-phiC); % deltaz = delta*sin(rem(x(4),2*pi)-phiC); deltay = delta*cos(x(4)-phiC); deltaz = delta*sin(x(4)-phiC); %--------------------------He thong luc --------------------------- hsl = 1/2*ro*S*V^2; hsm = hsl*l; matranquay = [1 0 0;0 cos(x(4)) -sin(x(4));0 sin(x(4)) cos(x(4));]; Lt = hsl*CLalpha*[0;-beta;-alpha;];% luc nang tong hop Dt = -hsl*CD*[1; beta; alpha;]; % luc can tong hop Kt = hsl*(Cypalpha_F*(x(11)*l/V)+Cypalpha_A*(x(12)*l/V))*[0;-alpha; beta;];% luc Magnus tong hop FC = hsl*CNdelta*(matranquay*[0;deltay;-deltaz;]-[0;beta;alpha;]); % luc dieu khien MC = hsm*CMdelta*(matranquay*[0;deltaz;deltay;]-[0;-alpha;beta;]); % mm dieu khien Mp = hsm*CMalpha*[0;alpha;-beta;]; % momen lat tong hop Mm = hsm*(-CMpalpha_F*x(11)*l/V-CMpalpha_A*x(12)*l/V)*[0;-beta;-alpha;];% mo men magnus tong hop Mpd = hsm*(CMq + CMalpha_dot)*l/V*[0;x(13);x(14);]; % momen giam chan tong hop Mrd_F = hsm*Clp_F*x(11)*l/V; % MM giam chan roll cho phan truoc Mrd_A = hsm*Clp_A*x(12)*l/V; % MM giam chan roll cho phan sau Ft = Lt + Dt + Kt + FC; Mt = MC + Mp + Mm + Mpd; X = Ft(1); Y = Ft(2); Z = Ft(3); Mxf = Mrd_F; Mxa = Mrd_A; My = Mt(2); Mz = Mt(3); %---------------He phuong trinh vp ------------------------------------------ dx(1) = X / m - g * sin(x(6)) + x(14) * x(2) - x(13) * x(3); dx(2) = Y / m - x(14) * x(1) - x(14) * tan(x(6)) * x(3); dx(3) = Z / m + g * cos(x(6)) + x(13) * x(1) + x(14) * tan(x(6)) * x(2); dx(4) = x(11) + x(14) * tan(x(6)); dx(5) = x(12) + x(14) * tan(x(6)); dx(6) = x(13); dx(7) = 0.1e1 / cos(x(6)) * x(14); dx(8) = cos(x(6)) * cos(x(7)) * x(1) - sin(x(7)) * x(2) + sin(x(6)) * cos(x(7)) * x(3); dx(9) = cos(x(6)) * sin(x(7)) * x(1) + cos(x(7)) * x(2) + sin(x(6)) * sin(x(7)) * x(3); dx(10) = -sin(x(6)) * x(1) + cos(x(6)) * x(3); dx(11) = 0.1e1 / Ixf * (Mxf + Cv * (x(12) - x(11))); P8 dx(12) = 0.1e1 / Ixa * (Mxa - Cv * (x(12) - x(11))); dx(13) = My / Iy - x(14) * (Ixa * x(12) + Ixf * x(11)) / Iy - x(14) ^ 2 * tan(x(6)); dx(14) = Mz / Iy + x(13) * (Ixa * x(12) + Ixf * x(11)) / Iy + x(13) * x(14) * tan(x(6)); dx(15) = FC(1); dx(16) = FC(2); dx(17) = FC(3); dx(18) = MC(1); dx(19) = MC(2); dx(20) = MC(3); end 3.3. Hệ số khí động đạn pháo 155mm function CM155 = hskddual155(M) A=[ 2.3271e-01 -1.0414e+00 1.4376e+00 -5.1653e-01 -1.1341e-01 2.4305e-01 5.6646e-02 4.9348e-02 -4.5156e+00 1.8257e+01 -2.1660e+01 9.0535e+00 1.1221e+01 -4.5882e+01 5.5675e+01 -1.3863e+01 -5.2544e+01 2.2223e+02 -2.8582e+02 8.9803e+01 -1.5522e-01 6.8767e-01 -9.6253e-01 2.9428e-01 -4.2739e+00 1.7781e+01 -2.2305e+01 7.3080e+00 -1.1612e-02 4.8302e-02 -6.0606e-02 1.9900e-02 -1.1622e+00 5.4148e+00 -7.8892e+00 2.4624e+00 -6.8322e-01 3.1579e+00 -4.5710e+00 1.4292e+00 5.4848e-01 -2.4597e+00 3.4064e+00 -8.8261e-01 4.8053e-01 -2.1499e+00 2.9774e+00 -8.1976e-01 4.6006e+00 -1.9140e+01 2.3878e+01 -7.5192e+00 1.2360e+01 -5.1267e+01 6.3736e+01 -2.0043e+01]; CM155 = A*[M^3; M^2; M; 1;]; end PHỤ LỤC 4 Chương trình tính toán thông số chuyển động bay của đạn phản lực 122mm lắp khoang điều khiển 4.1. Chương trình chính clc; clear all; close all; u0 = 42.7095; v0 = 0; w0 = 0; phif0 = 0; % goc quay ban dau cua phan phia truoc phia0 = 0; % goc quay ban dau cua phan phia sau theta0 = 45*pi/180; % goc ban ban dau psi0 = 0; %goc huong ban dau xe0 = 0; ye0 = 0; ze0 = 0; pf0 = 0;% rad/s toc do ban dau phan trc pa0 = 18.3342; % rad/s toc do ban dau phan sau q0 = 0; % toc do goc pitch ban dau r0 = 0; %toc do goc yaw ban dau x0 = [u0 v0 w0 phif0 phia0 theta0 psi0 xe0 ye0 ze0 pf0 pa0 q0 r0]; %index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 deltat = 0.01; global phiC deltaC tC dCNdelta dCMdelta hsCD hsCMalpha hsdamp deltamf hsCD = 0; % he luc can hsCMalpha = 0; %he so mo men chuc ngoc; Khao sat DeltaC da dieu chinh he so P9 len 2.5 den tang de theta<90do trong thoi gian khao sat 50s %Thay doi dai tu -60% den 20% dCMdelta = 0;%0.3;%0.6;%0.5;%0.3;%0.45;%0.5; % so gia dieu chinh CMdelta theo % %Thay doi khoi luong deltamf = 0; %kg - khoi khuong tang phan truoc dCNdelta = 0;%-0.15;%0.3;%0.6;%0.5;%0.3;%0.45;%0.5;% so gia dieu chinh CNdelta hsdamp = 0.5; % hs giam giam chan phiC = pi; % Pha ban dau cua dao dong canh lai, dv radian % phiC0 =0; % phiC180 = pi; % 0 --lech phai (duong OY) % pi --lech trai (am OY) % pi/2 --mui keo xuong duoi % 3pi/2 --mui dan keo len delta0 = 10;%5;%18;%12;%8;%12; % Bien do canh lai dv do deltaC = delta0*pi/180; % Bien do canh lai tC = 35; %thoi gian canh hoat dong Tf = 62; [T,X,Ts_xuat] = rungekutta(Tf,x0,deltat,@hptdualspinkdk21_ksIxf); u = X(:,1); v = X(:,2); w = X(:,3); phif = X(:,4); phia = X(:,5); theta = X(:,6); psi = X(:,7); xe = X(:,8); ye = X(:,9); ze = X(:,10); pf = X(:,11); pa = X(:,12); q = X(:,13); r = X(:,14); V = sqrt(u.^2+v.^2+w.^2); alpha = w./V; beta = v./V; alphat = sqrt(alpha.^2+beta.^2); ze_0= find(-ze>0); Tmax = max(ze_0)*deltat tam = xe(max(ze_0)) figure(1); subplot(2,2,1); plot(T,alpha*180/pi);title('alpha-goc tan'); grid on; subplot(2,2,2); plot(T,beta*180/pi);title('beta-goc truot canh');grid on; subplot(2,2,3); plot(T,V);title('V tong');grid on; subplot(2,2,4); plot(T,theta*180/pi);title('theta');grid on; figure(4); subplot(2,2,1); plot(T,60*pf/2/pi);title('pf - Toc do goc phan truoc [rps]'); subplot(2,2,2); plot(T,60*pa/2/pi);title('pa - Toc do goc phan sau [rps]'); subplot(2,2,3); plot(T,q);title('q - Toc do goc pitch [rad/s]'); subplot(2,2,4); plot(T,r);title('r - Toc do goc yaw [rad/s]'); figure(2); subplot(2,2,1); plot(T,-ze); title('ze theo T'); grid on; subplot(2,2,2); plot3(xe,ye,-ze); xlabel('Tam');ylabel('Huong');zlabel('Cao do'); subplot(2,2,3); plot(xe,-ze);title('ze theo xe'); subplot(2,2,4); plot(xe,-ye);title('ye theo xe'); figure(5) plot(T,180/pi*sqrt(alpha.^2+beta.^2)) figure(6) plot(180/pi*beta,180/pi*alpha) P10 van_toc_tinh_SM = [min(V(tC/deltat:round(Tmax)/deltat)) max(V(tC/deltat:round(Tmax)/deltat))]; figure(7); plot(T,180/pi*psi); 4.2. Hệ phương trình vi phân chuyển động bay đạn phản lực 122mm lắp khoang điều khiển function [dx,g] = hptdualspinkdk21_ksIxf(t,x) dx = zeros(14,1); global deltaC tC dCNdelta dCMdelta hsCD hsCMalpha hsdamp phiC deltamf hskl = 1; mf = hskl*10+deltamf; ma = (-11.4*t + 66.6)*(t=1.74); m = mf + ma; Ixf = hskl*0.02837; Ixa = (-0.0144*t+0.15)*(t=1.8); Iy = hskl*(-8.474*t+62.5428 )*(t=1.8); d = 0.122; l = d; S = 1/4*pi*d^2; g = 9.8; Cv = 0.01; % gia tri truoc 0.02 CR = 0.001; % Maro = khiquyen(-x(10)); %toc do am thanh va mat do theo do cao ro = 1.2; %mat do khong khi at = 330; %van toc am thanh % ro = Maro(1);% mat do am thanh theo do cao; % at = Maro(2); %vt am thanh theo do cao V = sqrt(x(1)^2+x(2)^2+x(3)^2); M = V/at; %hs Mach alpha = x(3)/x(1); beta = x(2)/V; %----------Hs khi dong ------------------------------------------------------ Adnm = hskdkdkbm21_ksIxf(M); CD = Adnm(1); CLalpha = Adnm(2); Cl_A = Adnm(4); Clp_F = Adnm(5); Clp_A = Adnm(6); CMalpha = Adnm(7)+hsCMalpha*Adnm(7); CMq_alpha_dot = hsdamp*Adnm(8); % co he so dieu chinh giam chan CNdelta = 2*Adnm(9)+2*(dCNdelta)*Adnm(9); % co dieu chinh them so gia, da tang them 30% de dam bao Do du tru on dinh tinh duoi 10% (tuong duong tang dien tich canh lai len 30%) CMdelta = 2*Adnm(10)+2*(dCMdelta)*Adnm(10); % co dieu chinh them so gia, da tang them 30% de dam bao Do du tru on dinh tinh duoi 10%(tuong duong tang dien tich canh lai len 30%) Cypalpha_F = 0; CMpalpha_F = 0; %---------Goc lat canh lai ------------------------------------------------- % delta = deltaC*(t>tC&&t<tCf); delta = deltaC*(t>tC); deltay = delta*cos(x(4)-phiC); deltaz = delta*sin(x(4)-phiC); %--------------------------He thong luc ------------------------------------- hsl = 1/2*ro*S*V^2; hsm = hsl*l; matranquay = [1 0 0;0 cos(x(4)) -sin(x(4));0 sin(x(4)) cos(x(4));]; P11 Lt = hsl*CLalpha*[0;-beta;-alpha;];% luc nang tong hop Dt = -hsl*CD*[1; beta; alpha;]; % luc can tong hop Kt = hsl*Cypalpha_F*(x(11)*l/V)*[0;-alpha; beta;];% luc Magnus chi cho phan ngoi F Fdc = 23000*(t<1.74)*[1;0;0;];% luc day dong co FC = hsl*CNdelta*(matranquay*[0;deltay;-deltaz;]-[0;beta;alpha;]); % luc dieu khien MC = hsm*CMdelta*(matranquay*[0;deltaz;deltay;]-[0;-alpha;beta;]); % mm dieu khien Mp = hsm*CMalpha*[0;alpha;-beta;]; % momen lat tong hop Mm = hsm*(-CMpalpha_F*x(11))*l/V*[0;-beta;-alpha;];% mo men magnus chi cho phan ngoi F Mpd = hsm*(CMq_alpha_dot)*l/V*[0;x(13);x(14);]; % momen giam chan tong hop Mrd_F = hsm*Clp_F*x(11)*l/V; % MM giam chan roll cho phan truoc Mrd_A = hsm*(Cl_A+Clp_A*x(12)*l/V); % MM giam chan roll cho phan sau Ft = Lt + Dt + Kt + Fdc + FC; Mt = MC + Mp + Mm + Mpd; X = Ft(1); Y = Ft(2); Z = Ft(3); Mxf = Mrd_F; Mxa = Mrd_A; My = Mt(2); Mz = Mt(3); %---------------He phuong trinh vp --------------------------- dx(1) = X / m - g * sin(x(6)) + x(14) * x(2) - x(13) * x(3); dx(2) = Y / m - x(14) * x(1) - x(14) * tan(x(6)) * x(3); dx(3) = Z / m + g * cos(x(6)) + x(13) * x(1) + x(14) * tan(x(6)) * x(2); dx(4) = x(11) + x(14) * tan(x(6)); dx(5) = x(12) + x(14) * tan(x(6)); dx(6) = x(13); dx(7) = 0.1e1 / cos(x(6)) * x(14); dx(8) = cos(x(6)) * cos(x(7)) * x(1) - sin(x(7)) * x(2) + sin(x(6)) * cos(x(7)) * x(3); dx(9) = cos(x(6)) * sin(x(7)) * x(1) + cos(x(7)) * x(2) + sin(x(6)) * sin(x(7)) * x(3); dx(10) = -sin(x(6)) * x(1) + cos(x(6)) * x(3); dx(11) = 0.1e1 / Ixf * (Mxf + Cv * (x(12) - x(11))); dx(12) = 0.1e1 / Ixa * (Mxa - Cv * (x(12) - x(11))); dx(13) = My / Iy - x(14) * (Ixa * x(12) + Ixf * x(11)) / Iy - x(14) ^ 2 * tan(x(6)); dx(14) = Mz / Iy + x(13) * (Ixa * x(12) + Ixf * x(11)) / Iy + x(13) * x(14) * tan(x(6)); end 4.3. Hệ số khí động đạn phản lực 122mm lắp khoang điều khiển function Adnm = hskdkdkbm21_ksIxf(M) Adnm=zeros(1,10); % CD CLalpha Cl_F Cl_A Clp_F Clp_A cmalpha cmalphaq+cmalpha_dot CNdelta CMdelta a = [0.2 0.53 20.5 0 0.101 -5.9 -8.84 -128.37 -2147.45 1.9478 17.7509 0.4 0.49 20.77 0 0.102 -6.12 -9.06 -129.85 -2202.25 2.0340 18.5381 0.6 0.46 21.24 0 0.105 -6.5 -9.35 -132.22 -2333.15 2.2341 22.3628 0.8 0.45 21.72 0 0.105 -7.12 -9.38 -131.08 -2545.6 2.7336 24.9174 1 0.54 23.76 0 0.108 -8.32 -9.5 -140.37 -3034.3 4.3255 35.6732 1.1 0.6 27.22 0 0.145 -9.83 -10.5 -170.91 -3699.8 P12 3.3774 30.6197 1.2 0.56 27.24 0 0.163 -8.87 -11.06 -184.93 -3978.9 2.4658 22.3497 1.3 0.6 24.66 0 0.178 -7.82 -11.5 -167.14 -4232.25 2.0126 18.2403 1.4 0.57 23.39 0 0.184 -6.48 -11.45 -168.44 -4204.35 1.6942 15.3540 1.5 0.55 21.39 0 0.18 -5.71 -10.33 -150.23 -3698.05 1.4693 13.3147 1.6 0.58 20 0 0.164 -5.21 -9.41 -138.77 -3785.7 1.3184 11.9466 1.8 0.52 16.99 0 0.15 -4.48 -8.58 -106.46 -2975.05 1.1067 10.0284 2 0.5 15.44 0 0.137 -3.95 -7.83 -91.49 -2752.7 0.9610 8.7075 2.2 0.48 14.21 0 0.126 -3.54 -7.21 -79.90 -2634.15 0.8524 7.7233]; if M<a(1,1) for j = 2:11 m = a(1,j)/a(1,1); Adnm(j-1)=m*M; end end if M>a(14,1) Adnm = a(14,2:11); end if M<=a(14,1) for i = 1:14 if M==a(i,1) Adnm = a(i,2:11); elseif M>a(i,1)&&M<a(i+1,1) for j = 2:11 k = (a(i+1,j)-a(i,j))/(a(i+1,1)-a(i,1)); Adnm(j-1)=k*(M-a(i,1))+a(i,j); end end end end PHỤ LỤC 5 Chương trình tính toán thông số ổn định đạn phản lực 122mm lắp khoang điều khiển clc; global deltaC phiC tC dCNdelta dCMdelta hsCMalpha hsdamp deltamf nT = length(T); Mt = V/330; d = 0.122; l = d; S = 1/4*pi*d^2; ro = 1.2; g = 9.8; Ltl = 3.32; % chieu dai ten lua %--------------------------------------------- %Xac dinh khoi luong m = zeros(nT,1); hskl = 1; % he so khao sat thay doi khoi luong Ixa = zeros(nT,1); Ixf = hskl*0.02837; Iy = zeros(nT,1); mf = hskl*10+deltamf; for j = 0:(nT-1) P13 t = deltat*j; i = j+1; m(i) = mf + (-11.4*t + 66.6)*(t=1.8); Ixa(i) = (-0.0144*t+0.15)*(t=1.8); Iy(i) = hskl*(-8.474*t+62.5428 )*(t=1.8); end %--------------------------------------------- %Xac dinh he so khi dong CD = zeros(nT,1); CLalpha = zeros(nT,1); Cl_A = zeros(nT,1); Clp_F = zeros(nT,1); Clp_A = zeros(nT,1); CMalpha = zeros(nT,1); CMq_alpha_dot = zeros(nT,1); CNdelta = 2*zeros(nT,1); CMdelta = 2*zeros(nT,1); Cypalpha_F = zeros(nT,1); CMpalpha_F = zeros(nT,1); for jj=1:nT temp_hskd=hskdkdkbm21_ksIxf(Mt(jj)); CD(jj) = temp_hskd(1); CLalpha(jj) = temp_hskd(2); Cl_A(jj) = temp_hskd(4); Clp_F(jj) = temp_hskd(5); Clp_A(jj) = temp_hskd(6); CMalpha(jj) = temp_hskd(7)+hsCMalpha*temp_hskd(7); CMq_alpha_dot(jj) = hsdamp*temp_hskd(8); CNdelta(jj) = 2*temp_hskd(9)+2*(dCNdelta)*temp_hskd(9); %da tang them 30% de dam bao Do du tru on dinh tinh duoi 10% (tuong duong tang dien tich canh lai len 30%) CMdelta(jj) = 2*temp_hskd(10)+2*(dCMdelta)*temp_hskd(10); %da tang them 30% de dam bao Do du tru on dinh tinh duoi 10% (tuong duong tang dien tich canh lai len 30%) Cypalpha_F(jj) = 0; CMpalpha_F(jj) = 0; end %------------------------------------ %Xac dinh cac he so gan sao----- hsS = ro*S*l/2/m(nT); % he so danh sao m lay khoi luong khi dc chay het CD_S = hsS*CD; CLalpha_S = hsS*CLalpha; Cl_A_S = hsS*Cl_A; Clp_F_S = hsS*Clp_F; Clp_A_S = hsS*Clp_A; CMalpha_S = hsS*CMalpha; CMq_alpha_dot_S = hsS*CMq_alpha_dot; CNdelta_S = hsS*CNdelta; CMdelta_S = hsS*CMdelta; Cypalpha_F_S = hsS*Cypalpha_F; CMpalpha_F_S = hsS*CMpalpha_F; %------------------- Fdc =zeros(nT,1);%luc day dong co for j = 0:(nT-1) t = deltat*j; jj = j+1; Fdc(jj) = 23000*(t<1.8); % luc day dong co end gS = g*l*sin(theta)./(V.^2); fdcS = Fdc*l/m(nT)./(V.^2); gS = gS-fdcS; % dieu chinh lai g* %------------------------------------------- % Xac dinh cac thong so trong dieu kien on dinh P14 kt = sqrt(Iy(nT)/(m(nT)*l*l)); kt_2 = kt^-2;% ban kinh hoi chuyen ngang epsilon = 0.5*tan(theta); deltaCt = zeros(nT,1); if tC <=Tf ntC = tC/deltat; deltaCt(ntC:nT) = deltaC*ones(nT-ntC+1,1); end %----------------------------------------- H = CLalpha_S - CD_S - kt_2*CMq_alpha_dot_S; M = kt_2*(CMalpha_S + CMdelta_S); P = l*(Ixf*pf+Ixa.*pa)./(Iy.*V); TT = CLalpha_S - gS +fdcS+ 1./P*kt_2.*CMpalpha_F_S.*pf*l./V; G = -1i*(CD_S+gS-fdcS+kt_2*CMq_alpha_dot_S+1i*P)*g*l.*cos(theta)./V.^2; % dang chia thu cho 2 cua CMq C = -kt_2*CMdelta_S.*deltaCt*exp(1i*phiC); zetae = -G./(M+1i*P.*TT)-C./(M+1i*P.*TT); alphae = imag(zetae); betae = real(zetae); ea = epsilon.*alphae; eb = epsilon.*betae; PT = P.*TT; H_S = (1-ea).*(H-2*gS); M_S = M./(1-ea) + eb.*PT./(1-ea).^2 + eb.*P.*M./(1-ea)./H_S... +((ea.*PT+eb.*M).^2-2*eb.*PT.*M)./(1-2*ea)./H_S./H_S; DT =(ea.*PT+eb.*M).^2./((1-2*ea).*(H_S.^2)); p1 = 2*(H-2*gS); p2 = (H-2*gS).^2 + P.^2 -(2-2*ea).*M./(1-2*ea) - 2*eb.*PT./(1-2*ea); Sg_S = P.^2/4./M_S; Sd_S = 2*TT./H_S; VT = 1./Sg_S; VP = Sd_S.*(2-Sd_S); %------------------------------ % Xac dinh bien do cua canh % DELTA = ((H-2*gS).*P.*M - 2*P.*TT.*M).^2-4*(M.*(H-2*gS).^2-(H-2*gS- TT).*TT.*P.^2).*M.^2; DELTA = (P.^2-4*M).*(H-2*gS).^2.*M.^2; B1=-((H-2*gS).*P.*M-2*P.*TT.*M+sqrt(DELTA))/2./M.^2; B2=-((H-2*gS).*P.*M-2*P.*TT.*M-sqrt(DELTA))/2./M.^2; neta = epsilon.*kt_2.* CMdelta_S./M; tdk = 34;%thoi diem qua dinh quy dao CHI CO TAC DUNG VE ndk = tdk/deltat; deltamax_1 = B2(ndk:nT)./neta(ndk:nT)/cos(phiC); deltamax_2 = B1(ndk:nT)./neta(ndk:nT)/cos(phiC); deltamax_temp2 = min(deltamax_2)*180/pi deltamax_temp1 = min(deltamax_1)*180/pi Ssd = VT-VP; %ve trai cua dieu kien on dinh 2 Sss = H - 2*gS-fdcS; T_max_od= find(Ssd<0,1,'last')*deltat figure(8); plot(T(ndk:nT),deltamax_1*180/pi,T(ndk:nT),deltamax_2*180/pi) ylabel('(p1;p2)');xlabel('thoi gian,[s]'); hold on; axis([0 53 -0.0002 0.0015]); grid on; %---------Bieu dien dieu kien on dinh theo thoi gian ------------------ figure(11); plot(T,zeros(1,nT),T,Ssd,'Linewidth',1.5,'Color','black'); ylabel('Thong so on dinh Ssd');xlabel('thoi gian,[s]'); hold on; axis([0 Tf -5e5 2e5]);grid on; %-----bieu dien on dinh trong he toa do 1/Sg, Sd--------------------------- figure(12) P15 plot(T,zeros(1,nT),T,Sss,'Linewidth',1.5,'Color','black'); ylabel('Thong so on dinh Sss');xlabel('thoi gian,[s]'); hold on; axis([0 Tf -0.001 0.005]);grid on; %--------------Luc va mo men khi dong tac dong len dan--------- hsl = 1/2*ro*S*V.^2; hsm = hsl*l; Lt = [hsl';hsl';hsl';].*[CLalpha';CLalpha';CLalpha';].*[zeros(1,nT);-beta';- alpha';];% luc nang tong hop Dt = [-hsl';-hsl';-hsl';].*[CD';CD';CD';].*[ones(1,nT); beta'; alpha';]; % luc can tong hop FC = [hsl';hsl';hsl';].*[CNdelta';CNdelta';CNdelta';].*([zeros(1,nT);deltaCt'.*one s(1,nT)*cos(phiC);deltaCt'.*ones(1,nT)*sin(phiC);]- [zeros(1,nT);beta';alpha';]); Ft = Lt+Dt+FC; MC = [hsm';hsm';hsm';].*[CMdelta';CMdelta';CMdelta';].*([zeros(1,nT);- deltaCt'.*ones(1,nT)*sin(phiC);deltaCt'.*ones(1,nT)*cos(phiC);]- [zeros(1,nT);-alpha';beta';]); % mm dieu khien Mp = [hsm';hsm';hsm';].*[CMalpha';CMalpha';CMalpha';].*[zeros(1,nT);alpha';- beta';]; % momen lat tong hop Mpt = Mp+MC; Fz = Ft(3,:); Dx = Dt(1,:); My = Mpt(2,:); % On dinh tinh SMt = (CMalpha+CMdelta)./(CLalpha+CNdelta); % Do du tru on dinh tinh theo thoi gian ndelta = 1/2*ro*S*V.^2.*(CNdelta-CLalpha.*CMdelta./CMalpha)./(m*g); figure(13); ntc = tC/deltat; % chi so cac vector tai thoi diem dieu khien SM = SMt(ntc:nT); plot(T(ntc:nT),SM,'Linewidth',1.5,'Color','black'); title('Do du tru on dinh tinh SM'); grid on; xlabel('thoi gian,[s]'); ylabel('Do du tru on dinh tinh,[duong kinh dan]'); title('Do du tru on dinh tinh SM'); grid on; SMmin = min(SM) SMmax = max(SM) SMmin_phantram = min(abs(SM)*d/Ltl*100) SMmax_phantram = max(abs(SM)*d/Ltl*100) figure(18) plot(T,alphae*180/pi,T,betae*180/pi,'-- ','Linewidth',1.5,'Color','black');legend('alphae','betae'); grid on; axis([21 Tf -1 8]) % plot(T,zeros(1,length(T)),T,SSD0,T,SSD1,T,SSD1_5,T,SSD2,'Linewidth',1.5,'Colo r','black'); ylabel('Goc tan va goc truot canh can bang,[do]'); xlabel('thoi gian,[s]'); %------------------------------------ alphae_do = alphae*180/pi; alphae_do_tb = mean(alphae_do(ntc:nT)) %goc tan can bang trung binh CE = alphae_do_tb/delta0 % chi so hieu qua dieu khien %------------- Xac dinh L/D--------------------------- Le = [hsl';hsl';hsl';].*[CLalpha';CLalpha';CLalpha';].*[zeros(1,nT);-betae';- alphae';];% luc nang tong hop De = [-hsl';-hsl';-hsl';].*[CD';CD';CD';].*[ones(1,nT); betae'; alphae';]; % luc can tong hop Let = zeros(nT,1); Det = zeros(nT,1); LDe = zeros(nT,1); for ii =1:nT Let(ii) = sqrt(Le(1,ii)^2+Le(2,ii)^2+Le(3,ii)^2); P16 Det(ii) = sqrt(De(1,ii)^2+De(2,ii)^2+De(3,ii)^2); LDe(ii) = Let(ii)/Det(ii); end figure(20); plot(T,CD_S); figure(21); plot(T,gS);