Về một phương pháp xây dựng mô hình thuật toán phương pháp phần tử hữu hạn tính toán ứng suấtbiến dạng thân vỏ tên lửa đối hạm KH-35E

Tên lửa & Thiết bị bay N.T. Bình, N.M. Tuấn, P.T. Lai, “Về một phương pháp xây dựng mô hình KH35-E.” 26 VỀ MỘT PHƯƠNG PHÁP XÂY DỰNG MÔ HÌNH THUẬT TOÁN PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TÍNH TOÁN ỨNG SUẤT- BIẾN DẠNG THÂN VỎ TÊN LỬA ĐỐI HẠM KH-35E Nguyễn Thanh Bình1*, Nguyễn Minh Tuấn2, Phan Tương Lai 2 Tóm tắt: Bài báo sử dụng lý thuyết tấm mỏng, phương pháp phần tử hữu hạn để nghiên cứu phát triển và xây dựng thuật toán xác định trường ứng suất - biến dạng của thân vỏ tên lửa đối hạm Kh-35E. Nghiên cứu thiết lập mô hình bài toán đối với thân vỏ tên lửa được rời rạc hóa bằng các phần tử tấm phẳng dạng tam giác liên tục chịu tải trọng khí động, xây dựng thuật toán đối với phần tử tam giác 3 nút, phân tích lựa chọn hàm xấp xỉ chuyển vị, xây dựng các biểu thức biểu diễn các quan hệ chuyển vị, quan hệ biến dạng, quan hệ ứng suất theo chuyển vị và dựa vào phiến hàm thế năng toàn phần tối thiểu, xác định ma trận độ cứng phần tử, từ đó thiết lập trạng thái ứng suất tại tâm phần tử. Từ khóa: Tên lửa đối hạm, Phương pháp phần tử hữu hạn, Lý thuyết tấm và vỏ. 1. ĐẶT VẤN ĐỀ Kết cấu thân cánh tên lửa là một hệ thống thống nhất vừa tạo lực nâng vừa chịu tương tác với môi trường khí quyển khi bay và làm việc trong điều kiện hết sức phức tạp và ngặt nghèo. Tính toán độ bền kết cấu thân cánh là một trong những nhiệm vụ quan trọng hàng đầu trong quá trình thiết kế, chế tạo, thử nghiệm tên lửa nói chung và tên lửa hành trình đối hạm nói riêng. Việc nghiên cứu riêng lẻ từng quá trình bằng sự ứng dụng sự phát triển vượt bậc của công nghệ thông tin và các mô hình toán - lý có thể khảo sát kỹ càng các mô hình đối tượng có hình dạng phức tạp kết hợp với ảnh hưởng đồng thời của rất nhiều yếu tố bằng phương pháp tính toán, mô phỏng số trên máy tính. Khi ứng dụng phương pháp phần tử hữu hạn, bề mặt cong liên tục của thân tên lửa được tập hợp bởi các phần tử tấm tam giác nhỏ, dẫn đến giải bài toán ứng suất phẳng và uốn đối với thân tên lửa. Xấp xỉ bề mặt thân và cánh tên lửa Kh35-E trong quá trình bay hành trình và tăng tốc trong giai đoạn cuối bằng các phần tử tấm phẳng dạng tam giác được chỉ ra trên hình 1. Khung thân Khung thân 1 2 3 w y x w* y* x* 1 2 3 w y x w* z* x* y* z a) b) Phân tu tam giác , Hình 1. Xấp xỉ bề mặt thân và cánh tên lửa Kh35-E bằng tổ hợp phần tử tấm phẳng dạng tam giác. 2. XÂY DỰNG THUẬT TOÁN PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN XÁC ĐỊNH TRƯỜNG ỨNG SUẤT-BIẾN DẠNG THÂN TÊN LỬA KH35-E Các giả thiết: - Chiều dày thân là hằng số và được coi như là vỏ mỏng. - Vỏ bọc của thân tên lửa được gắn chặt với khung thân ghép nối các khoang và thanh dọc tăng cường của thân tên lửa. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 27 - Bỏ qua ứng suất theo hướng z, tức là z = 0 - Không xét đến lực cắt, tức là xy = yz =0 Với những giả thiết trên, mỗi phần tử hữu hạn trong kết cấu vỏ chịu các thành phần tải trọng tác dụng độc lập: lực màng và lực vuông góc với mặt phẳng vỏ. Xét một phần tử điển hình trong kết cấu, mỗi nút phần tử chịu hai chuyển vị màng u và v, độ võng w, góc xoay x , y , góc xoắn z trong mặt phẳng. Lực tác dụng lên phần tử điển hình bao gồm hai lực màng Fx , Fy , hai mômen uốn Tx , Ty , lực pháp tuyến Fz , mômen xoắn Tz tương ứng với góc xoay z . Nhận thấy rằng, góc xoắn trong mặt phẳng z là rất nhỏ và có thể bỏ qua, đồng thời không ảnh hưởng đến các số hạng khác trong ma trận độ cứng phần tử. Toàn bộ chuyển vị và tải trọng tại nút 1 có thể được viết dưới dạng véctơ sau đây, [7]:                   1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 và T Tm u z x y z T Tm u z x y z x y z u v F F F T F F F T T T           (1) Trong đó, các chỉ số trên m và u ký hiệu trạng thái màng và trạng thái uốn tương ứng. Xác định các ma trận độ cứng phần tử và các phương trình cơ bản giải bài toán tính trường ứng suất - biến dạng thân tên lửa được xấp xỉ bằng tập hợp các phần tử phẳng tam giác được trình bày dưới đây. 2.1. Xác định ma trận độ cứng phần tử tam giác đối với trạng thái màng 1. Chọn hàm chuyển vị [f(x,y)] và xác định véctơ chuyển vị  ( , )mx y tại điểm bất kỳ trong mặt phẳng: 1 2 3 4 5 6 u x y v x y              (2) và hàm chuyển vị có thể được viết như sau, [8]:      ( , ) ( , ) 1 2 3 4 5 6 1 0 0 0 0 0 0 1 Tm m x y x y u x y f v x y                         (3) 2. Biểu diễn chuyển vị  ( , )mx y tại điểm bất kỳ bên trong phần tử theo chuyển vị nút  me Đặt các giá trị của tọa độ nút ( 1, 1) ( 3, 3),..., m m x y x yf       vào phương trình (3) sau đó giải để tính véctơ {} (véctơ hệ số chưa biết của đa thức được cho trong phương trình (2)):           11 1 21 1 1 32 2 2 2 2 4 3 3 3 5 3 3 6 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 m me m m m x y x y x y A x y x y x y                                                              (4) Như vậy, phương trình (4) xác định ma trận [Am] cho trường hợp cụ thể của phần tử tam giác đối với trạng thái màng. Từ phương trình (4) nhận được: {} = [Am]-1 {me} Phương trình tổng quát (3) được viết như sau:     1 ( , ) ( , ) m m m me x y x yf A           (5) 3. Biểu diễn quan hệ biến dạng  ( , )mx y tại điểm bất kỳ theo chuyển vị  ( , )mx y và theo chuyển vị nút  me , véctơ biến dạng được viết như sau, [8]:    ( , ) Tm x y x y xy    (6) Tên lửa & Thiết bị bay N.T. Bình, N.M. Tuấn, P.T. Lai, “Về một phương pháp xây dựng mô hình KH35-E.” 28 Từ quan hệ biến dạng - chuyển vị và đặt u, v từ phương trình (2) và lấy đạo hàm theo các biến x, y, nhận được:      ( , ) 2 6 3 5 T Tm x y x y xy          (7) hoặc có thể biểu diễn như sau:         1 ( , ) m m m m me m me x y C C A B                      (8) ma trận biến dạng màng: 2 3 3 1 1 2 3 2 1 3 2 1 3 2 2 3 1 3 3 1 2 1 1 2 1 0 0 0 1 0 0 0 2 m m m y y y y y y x x x x x x x x y y x x y y x x y y B C A                                  4. Quan hệ ứng suất  ( , )mx y theo biến dạng  ( , )mx y và chuyển vị nút  me :       ( , ) ( , ) Tm m x y x y xy x yA      hay      ( , )mx y m meA B       (9) trong đó, [A] - ma trận đàn hồi của trạng thái màng 5. Biểu diễn quan hệ ứng suất  ( , )mx y với tải trọng nút tĩnh tương đương, quan hệ lực nút với chuyển vị nút và nhận được ma trận độ cứng phần tử [Kme]. Ma trận độ cứng [Kme] là ma trận đối xứng với các thành phần ij ij mk k 11 21 31 41 51 61 21 22 32 42 52 62 31 32 33 43 53 63 41 42 43 44 54 64 51 52 53 54 55 65 61 62 63 64 65 66 1 4                         m m m m m m m m m m m m m m m m m m me m m m m m m m m m m m m m m m m m m k k k k k k k k k k k k k k k k k k K k k k k k k k k k k k k k k k k k k (10) 6. Thành lập ma trận chuyển vị ứng suất [Hme] Quan hệ ứng suất - chuyển vị:   ( , )m m mex y H     (11) Trong đó: [Hm] = [Am][Bm] Từ (11) có thể xác định được ứng suất  ( , )mx y tại điểm bất kỳ (x,y) trong phần tử. Các ứng suất nhận được có chứa các số hạng phụ thuộc theo tọa độ x và y cho nên để nhận được ứng suất tại một điểm nào đó trong phần tử các tọa độ của điểm đó phải được đặt vào trong ma trận [Hm]. Trong trường hợp phần tử tam giác lấy tâm tam giác làm tọa độ của x, y tức là ứng suất nhận được tại tâm của tam giác. 2.2. Xác định ma trận độ cứng của phần tử tam giác phẳng dạng tam giác chịu uốn [Kue] Tại mỗi nút của phần tử có một độ võng w và 2 góc xoay x và y theo các trục x và y tương ứng được xem là các bậc tự do của mỗi nút. Các góc xoay này chính là đạo hàm của hàm độ võng w theo y và x. Như vậy, 9 thành phần chuyển vị phải được xét cho mỗi phần tử có 9 bậc tự do, được viết dưới dạng véctơ sau đây, [7]:    1 1 1 2 2 2 3 3 3w w w Tue x y x y x y       (12) 1. Chọn hàm chuyển vị w(x,y) và xác định véctơ chuyển vị  ( , )x y tại điểm bất kỳ của phần tử, hàm độ võng w(x,y) được xấp xỉ bằng một đa thức chứa 9 tham số có dạng [7]: 2 2 3 2 2 3 1 2 3 4 5 6 7 8 9w( , ) ( )x y x y x xy y x x y xy y                  (13) hay được viết dưới dạng ma trận sau Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 29     2 2 3 2 2 3w( , ) 1 ( ) ( , )x y x y x xy y x x y xy y P x y      (14) trong đó,    1 2 3 4 5 6 7 8 9 T T           Đến đây, chúng ta kiểm tra tính tương thích của phần tử đang xét với hàm độ võng w(x,y) . Giả sử xét cạnh biên ij là cạnh biên chung giữa hai phần tử A và B kề sát nhau. Trong hệ tọa độ địa phương cạnh ij có phương trình y=0. Chuyển vị theo phương z hay độ võng theo cạnh này là: 2 3 ij 0 1 2 4 7w (w)y x x x        (15) Độ dốc theo phương x hay w x   dọc theo cạnh này là: 2 2 4 7 ij 0 w w 2 3 y x x x x                      (16) Vì nút i và j là nút chung của cả hai phần tử A và B nên các bậc tự do hay các chuyển vị nút của 2 nút này cũng là chung đối với hai phần tử có cạnh biên chung ij này. Tuy nhiên, nếu để ý tới 4 chuyển vị nút 0 i j w w w , , w , lim                i j xx x  ta có thể thấy rằng bằng phương pháp đồng nhất 4 bậc tự do này với giá trị hàm w và giá trị đạo hàm w         x tại hai điểm nút i (x=0) và j (x=a), tức là chúng có 4 điều kiện sau: wi= wj= wij ; w w w i j ij x x x                                tại x=0 và x=a. Từ đó, hoàn toàn có thể xác định được 4 tham số 1 , 2 , 4 và 7 một cách duy nhất. Do đó, w và w x   hoàn toàn xác định trên cạnh biên ij và tính liên tục của chuyển vị và độ dốc dọc theo cạnh này được bảo đảm khi chuyển từ phần tử A sang phần tử B. Độ dốc theo phương y hay w y   dọc theo cạnh biên chung ij (y =0): 23 5 8 ij 0 w w y y y x x                       (17) Có thể thấy rằng để xác định w y   dọc theo cạnh biên chung ij, ta cần xác định được 3 tham số 3 , 5 và 8 . Nhưng chúng ta chỉ có được hai phương trình từ việc thay thế vào hai điều kiện: w w y yi ij                 tại x = 0 và w w y yj ij                 tại x = a Do đó độ dốc vuông góc với biên ij là không xác định. Độ dốc w y   này tại các nút i và j là như nhau đối với cả hai phần tử, nhưng có thể là khác nhau tại các điểm khác dọc theo cạnh biên ij. Do vậy, phần tử này là không tương thích. Tuy nhiên, phần tử này vẫn có thể áp dụng được vì trên thực tế người ta thường sử dụng một loại phần tử tấm phẳng chỉ đòi hỏi sự liên tục của chuyển vị và góc xoay tại các nút. Do đó các bậc tự do của mỗi nút sẽ là giá trị của w, w x và w y tại nút đó. Tên lửa & Thiết bị bay N.T. Bình, N.M. Tuấn, P.T. Lai, “Về một phương pháp xây dựng mô hình KH35-E.” 30 2. Biểu diễn chuyển vị  ( , )ux y tại điểm bất kỳ bên trong phần tử theo chuyển vị nút  ue Đạo hàm (14) theo x và y nhận được véctơ chuyển vị, [7]:   2 2 3 2 2 3 2 2 2 2 1 ( ) 0 0 1 0 2 0 ( 2 ) 3 0 1 0 2 0 3 (2 ) 0 x y w x y x xy y x x y xy y x y x xy y x y x xy y                            (18) Thay các tọa độ nút: nút 1 (0,0); nút 2 (x2 , 0); nút 3 (x3 , y3) nhận được:    ue nA     (19) trong đó:  nA    - véctơ hệ số; nA   - ma trận hằng số Giải phương trình (19) nhận được:     1u ueA       (20) Đặt (20) vào (14) nhận được w = [P(x,y)] [Au]-1{ue} (21) 3. Biểu diễn  ( , )ux y tại điểm bất kì theo chuyển vị nút {ue} Các biến dạng trong tấm là , ,x y xy   được biểu diễn dưới dạng ma trận như sau:     2 2 2 ( , ) 2 2 w w w 2 T Tu x y x y xy z x y x y                   (22) Từ (21) và (22) véctơ biến dạng  ( , )ux y được biểu diễn qua chuyển vị nút {ue} như sau:           2 2 2 1 ( , ) 2 2 ( , ) ( , ) 2 ( , )               T u u ue x y z P x y P x y P x y A x yx y   (23) Hoặc có thể viết    ( , )u u uex y B     Trong đó, ma trận biến dạng uốn [Bu] là [Bu] = -z[Cu][Au]-1     2 2 2 2 2 0 0 0 2 0 0 6 2 0 ( , ) 0 0 0 0 0 2 0 2 6 0 0 0 0 2 0 0 (4 4 ) 0 2                             u x x y C P x y x y y x y x y (24) 4. Xác định ma trận độ cứng uốn của phần tử [Kue] Ma trận độ cứng uốn [Kue] phần tử tam giác 3 nút được xác định, [7]:      /2 1 12/2 TT t Tue u u u u u u t Ve S K B D B dV z dz A C D C A dS                                  (25) Vì [Au]-1 chỉ chứa các hằng số nên đưa ra khỏi dấu tích phân và nhận được     3 1 1 12 T Tue u u u u S t K A C D C dS A                            (26) Trong đó, t và S lần lượt là chiều dày và diện tích phần tử tam giác. Công thức (26) có thể viết như sau    1 1 T ue u uK A I A              (27) Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 31 Từ biểu thức (27) nhận được ma trận độ cứng của phần tử tam giác [Kue]. Ma trận [Kue] nhận được bằng cách nhân từng đôi các ma trận [I][Au]-1 sau đó nhân với ([Au]-1)T nhận được ma trận có kích thước (9x9) và ma trận có tính chất đối xứng ij  u u jik k . 5. Thành lập ma trận chuyển vị -ứng suất [Hue] Quan hệ ứng suất - chuyển vị được cho trong phương trình (23) và (24) như sau    ( , )u ue uex y H    Trong trường hợp phần tử tam giác, nếu lấy tâm tam giác làm tọa độ của x, y thì nhận được ứng suất tại tâm tam giác. 2.3. Xác định ma trận độ cứng tổng thể [Ke] của phần tử tam giác đồng thời chịu lực màng và chịu lực uốn Ma trận độ cứng phần tử [Ke] có kích thước (18x18) và được tập hợp bởi hai ma trận độ cứng của trạng thái màng và ma trận độ cứng của trạng thái uốn. Quan hệ giữa tất cả các lực và chuyển vị tại nút 1 như sau:             1 111 1 1 11 1 1 1 0 0 0 0 0 0 0                                    m mm u u u z z F k F F k T    (28) Đến đây, quan hệ giữa lực và chuyển vị tại nút 1 đã được xác định. Quan hệ tương tự như vậy sẽ tồn tại tại 2 nút còn lại của phần tử. Chẳng hạn như, lực tại nút 2 có thể biểu diễn theo chuyển vị tại nút 3 như sau:             2 323 2 2 23 3 2 3 0 0 0 0 0 0 0                                    m mm u u u z z F k F F k T    (29) Tập hợp tất cả các quan hệ giữa lực và chuyển vị nút cho tất cả các phần tử có thể viết dưới dạng phương trình (34) đối với bài toán tính toán trường ứng suất - biến dạng của thân tên lửa được xấp xỉ bằng những phần tử phẳng dạng tam giác và ma trận độ cứng [Ke] có kích thước (18x18) 11 21 31 11 21 31 21 22 32 21 22 32 31 32 33 31 32 33 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                                                                    m m m u m m m m m e u u m m m m u u u k k k k k k k k k K k k k k k k k k k 0 0 0 0 0 0 0 0 0 0                                    (30) Đối với những vị trí đặt khoang thân và thanh dọc của thân tên lửa thì ma trận độ cứng [Ke] trong (30) được tính như sau: với giả thiết vỏ bọc của thân tên lửa được gắn chặt với khung thân ghép nối và thanh dọc tăng cường của thân tên lửa, như vậy, với các vị trí liên kết này giả thiết các khung thân và thanh dọc được coi như dầm chịu lực màng và đồng Tên lửa & Thiết bị bay N.T. Bình, N.M. Tuấn, P.T. Lai, “Về một phương pháp xây dựng mô hình KH35-E.” 32 thời chịu uốn. Ma trận độ cứng phần tử đối với các vị trí khung thân và thanh dọc chịu lực được cộng thêm độ cứng đối với trạng thái màng và uốn của phần tử dầm, [7]: - Đối với trạng thái màng ij ij E         m Ak k L - Đối với trạng thái uốn ij ij 3 12EJ         uk k L - đối với vị trí độ võng w; ij ij 2 6EJ         uk k L - đối với vị trí góc xoay y ; Trong đó, A - Tiết diện mặt cắt ngang của khung thân hoặc thanh dọc; J - Mômen quán tính của mặt cắt ngang khung thân hoặc thanh dọc; L - Chiều dài phần tử dầm của khung thân hoặc thanh dọc. 2.4. Chuyển hệ trục tọa độ Khi xây dựng các công thức tính ma trân độ cứng của phần tử, thường chúng ta sử dụng một hệ tọa độ địa phương phù hợp sao cho việc xây dựng công thức là dễ dàng và đơn giản. Các hệ tọa độ địa phương này có thể là khác với hệ tọa độ tổng thể chung cho toàn hệ. Do đó, trước khi ghép nối phần tử để xây dựng ma trận độ cứng tổng thể, ta cần chuyển ma trận độ cứng trong hệ tọa độ địa phương vào hệ tọa độ tổng thể. Trong bài toán vỏ cũng như bài toán thanh, các phần tử kề nhau được nối nghiêng với nhau. Bởi vậy, trước khi tập hợp các phần tử và ma trận độ cứng của toàn hệ kết cấu được thành lập thì chuyển vị và lực của mỗi phần tử riêng biệt phải được biểu diễn trong hệ tọa độ tổng thể. Ma trận độ cứng phần tử đã nhận được trong hệ tọa độ tổng thể và có thể viết dưới dạng chuẩn đối với quan hệ lực và chuyển vị trong hệ tọa độ tổng thể như sau    * * *   e e eF K  (31) 3. MÔ PHỎNG TRƯỜNG ỨNG SUẤT - BIẾN DẠNG THÂN TÊN LỬA KH35-E SỬ DỤNG PHẦN MỀM ANSYS ANSYS sử dụng phương pháp phần tử hữu hạn để giải bài toán kết cấu thông qua việc rời rạc hóa các phần tử kết cấu thành dạng lưới và giải hệ phương trình (31) cho các phần tử lưới đó. Cấu trúc lưới trong ANSYS có nhiều dạng nhưng nghiên cứu sử dụng dạng lưới tứ diện để sát với kiểu xấp xỉ bề mặt thân tên lửa bằng các phần tử tấm phẳng dạng tam giác. Về cấu trúc không gian, lưới có cấu tạo kiểu tứ diện. Mô hình lưới thân tên lửa được chia trong ANSYS như hình 2. Do kết cấu và tải trọng tác dụng đối xứng nên tác giả tiến hành nghiên cứu cho ½ mô hình thân tên lửa. Phụ thuộc vào tài nguyên máy tính nên độ mịn của lưới chỉ đạt mức vừa phải, đáp ứng được độ hội tụ cần thiết của kết quả tính toán. Thiết lập các thuộc tính vật liệu cho các phần tử kết cấu là hợp kim nhôm AMg-6. Thiết lập kiểu liên kết giữa các bề mặt tiếp xúc dạng gắn chặt (bonded). Lựa chọn trường hợp tính toán khi tên lửa bay trong mặt phẳng thẳng đứng ở cuối giai đoạn phóng. Hình 2. Mô hình các phần tử lưới tứ diện thân tên lửa. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 33 Theo [1], vận tốc tên lửa đạt cực đại (ứng với số Mach ≈ 0.952). Tên lửa bay với góc tấn α = 00. Tải trọng tác dụng lên thân tên lửa gồm: - Lực đẩy động cơ tác dụng lên bề mặt trong của loa phụt T = 79480 (N). - Tải trọng khí động tác dụng lên bề mặt chảy bao. Từ bài báo [2] nhận được giá trị tải trọng khí động gồm lực cản và lực nâng như sau: Fx_kd = 5988 (N); Fy_kd = -780.8 (N). -Tải trọng gây ra do khối lượng (gồm trọng lực và lực quán tính), gọi là tải trọng khối lượng. Giả thiết khối lượng phân bố đều, tải trọng khối lượng tác dụng lên toàn bộ kết cấu thân tên lửa. Trong trường hợp tính, tải trọng khối lượng biểu diễn dưới dạng hai thành phần dọc trục theo Ox và vuông góc theo Oy tương ứng: * * ; * * ;x yS n m g P n m g    Trong đó: - xn , y n - hệ số quá tải theo các phương Ox và Oy: _ ; * x kd x T F n m g   _ ; * y kd y F n m g  - m = 592,3 (kg) – khối lượng tên lửa tại thời điểm xét; - g là gia tốc trọng trường. Thiết lập tải trọng tác dụng lên ½ thân tên lửa như hình 3. Giải bài toán nhận được trường phân bố ứng suất – biến dạng của thân tên lửa như hình 4 và hình 5. Hình 4. Trường phân bố ứng suất. Hình 5. Trường phân bố biến dạng. Nhận thấy ứng suất cực đại xuất hiện tại đáy khoang chứa động cơ phóng có giá trị σmax = 1,11x108 (Pa) nhỏ hơn nhiều so với ứng suất tới hạn của vật liệu hợp kim nhôm AMg-6 (σb = 3,2x10 8 Pa). Biến dạng tỉ đối cực đại εmax = 1,57x10 -3 cũng nhỏ hơn so với giá trị giới hạn là 2x10-3 nên kết cấu đảm bảo độ bền. 4. KẾT LUẬN Bài báo nghiên cứu thiết lập mô hình bài toán đối với thân vỏ tên lửa được rời rạc hóa bằng các phần tử tấm phẳng dạng tam giác liên tục chịu tải trọng khí động, đã xây dựng thuật toán đối với phần tử tam giác 3 nút, mỗi nút 6 bậc tự do và mỗi phần tử 18 bậc tự do, đã tiến hành phân tích lựa chọn hàm xấp xỉ chuyển vị, xây dựng các biểu thức biểu diễn các quan hệ chuyển vị, quan hệ biến dạng, quan hệ ứng suất theo chuyển vị và dựa vào phiến hàm thế năng toàn phần tối thiểu đã xác định được ma trận độ cứng phần tử vừa chịu lực màng và uốn đống thời, từ đó thiết lập được trạng thái ứng suất tại tâm phần tử. Hình 3. Tải trọng tác dụng lên ½ thân tên lửa. Tên lửa & Thiết bị bay N.T. Bình, N.M. Tuấn, P.T. Lai, “Về một phương pháp xây dựng mô hình KH35-E.” 34 Quá trình xây dựng thuật toán, các biểu thức quan hệ và phương trình được trình bày bằng ngôn ngữ véctơ và ma trận rõ ràng, thuận lợi cho việc lập trình trên máy tính. Bên cạnh đó bài báo cũng trình bày kết quả nghiên cứu sử dụng phần mềm ANSYS dựa trên mô hình toán đã nêu để xác định trường ứng suất - biến dạng thân tên lửa Kh-35E ở thời điểm cuối giai đoạn phóng. Thân tên lửa đảm bảo độ bền trong điều kiện tính toán. Việc phát triển xây dựng thuật toán phương pháp phần tử hữu hạn có ý nghĩa rất quan trọng, làm cơ sở vững chắc cho việc khai thác phần mềm Ansys để giải các bài toán đặt ra một cách chính xác, tin cậy. TÀI LIỆU THAM KHẢO [1]. Trần Đức Trung, Cao Bá Ninh, “Tính toán thiết kế mẫu nguyên lý và triển khai chế tạo một số cụm chức năng của tên lửa hành trình đối hải với vận tốc bay dưới âm”, (2005). [2]. T.M.Tuân, N.P.Thắng C.D.Lành, “Xác định một số thông số khí động thiết bị bay trong môi trường ANSYS CFX,” TC, Nghiên cứu KHCNQS, Đặc san CNTT (04- 2014), tr, 92-99. [3]. Nguyễn Hoa Thịnh, Hoàng Xuân Lượng, Nguyễn Đức Cương, Trịnh Hồng Anh, Nguyễn Minh Tuấn, “Kết cấu và tính toán độ bền khí cụ bay ”, Nhà xuất bản KHKT, (2005). [4]. Лизин В.Т., Пяткин В.А, “Проектирование тонкостенных конструкций. Машиностроение”, М., (1994). [5]. Оболенский Е.П., Сахаров Б.И., Сибиряков В.А. “Прочность летательных аппаратов и их агрегатов”. Машиностроение, М., (1995). [6]. Zienkiewics O.C., “The Finite Element Methods in Engineering Science”. Mc Graw- Hill, (1979). [7]. Petyt M., “Introduction to Finite Element Vibration Analysis”. Cambridge University Press, United Kingdom, (1998). [8]. Timosenko S.P., “Theory of Plane and Shells”. Mc Graw-Hill, (1951). ABSTRACT ALGORITHM OF SETTING UP MODEL FINITE ELEMENT METHOD TO CALCULATION STRESS - DEFORMATION OF KH-35E SHIP MISLES’S HOUSING This paper presents the method and the general scheme of the model identification of aerodynamics characteristics derivatives with respect to the values of the angular velocities and accelerations of the aircraft. To exclude the confusion we call six dimensionless transient characteristics of the aircraft relatively rigidly connected with the non-deformable part of the aircraft coordinate system. The development of iterative algorithms of identification of aerodynamic characteristics and the coefficients are based on the high-precision motion model aircraft, with the possibility of monitoring the results of calculations. Keywords: Ship misles, Finite Element Methods, Theory of Plane and Shells. Nhận bài ngày 09 tháng 9 năm 2015 Hoàn thiện ngày 01 tháng 12 năm 2015 Chấp nhận đăng ngày 25 tháng 12 năm 2015 Địa chỉ: 1Trung tâm CNCKCX - Viện KH-CNQS; 2Viện KH-CNQS; *Email: binhctm@yahoo.com