Tên lửa & Thiết bị bay<br />
<br />
VỀ MỘT PHƯƠNG PHÁP XÂY DỰNG MÔ HÌNH THUẬT TOÁN<br />
PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TÍNH TOÁN ỨNG SUẤT-<br />
BIẾN DẠNG THÂN VỎ TÊN LỬA ĐỐI HẠM KH-35E<br />
Nguyễn Thanh Bình1*, Nguyễn Minh Tuấn2, Phan Tương Lai2<br />
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 để<br />
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<br />
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<br />
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<br />
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,<br />
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<br />
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<br />
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<br />
trạng thái ứng suất tại tâm phần tử.<br />
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ỏ.<br />
<br />
1. ĐẶT VẤN ĐỀ<br />
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<br />
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<br />
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<br />
đầ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<br />
đố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<br />
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ô<br />
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<br />
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<br />
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<br />
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.<br />
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<br />
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.<br />
, Khung thân<br />
Khung thân Phân tu tam giác<br />
<br />
<br />
<br />
<br />
y* z* y* y 2<br />
y z<br />
2 x x<br />
3<br />
3 w1<br />
w<br />
1 x* x*<br />
* *<br />
w w<br />
a) b)<br />
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.<br />
2. XÂY DỰNG THUẬT TOÁN PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN XÁC ĐỊNH<br />
TRƯỜNG ỨNG SUẤT-BIẾN DẠNG THÂN TÊN LỬA KH35-E<br />
Các giả thiết:<br />
- Chiều dày thân là hằng số và được coi như là vỏ mỏng.<br />
- 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<br />
dọc tăng cường của thân tên lửa.<br />
<br />
<br />
26 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.”<br />
Nghiên cứu khoa học công nghệ<br />
<br />
- Bỏ qua ứng suất theo hướng z, tức là z = 0<br />
- Không xét đến lực cắt, tức là xy = yz =0<br />
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<br />
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ử<br />
đ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<br />
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<br />
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<br />
ứ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ỏ<br />
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ử.<br />
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]:<br />
T T<br />
1 1m 1u z1 u1 v1 x1 y1 z1<br />
T<br />
(1)<br />
T<br />
và F1 F 1<br />
m<br />
F 1<br />
u<br />
Tz1 F<br />
x1 Fy1 Fz1 Tx1 Ty1 Tz1<br />
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.<br />
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<br />
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<br />
giác được trình bày dưới đây.<br />
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<br />
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<br />
u 1 2 x 3 y<br />
mặt phẳng: (2) và hàm chuyển vị có thể được viết như sau, [8]:<br />
v 4 5 x 6 y<br />
u 1 x y 0 0 0 T<br />
f<br />
m<br />
( x, y )<br />
m<br />
( x, y)<br />
<br />
v 0 0 0 1 x y 1 2 3 4 5 6 (3)<br />
<br />
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 <br />
Đặt các giá trị của tọa độ nút (mx1, y1) ,..., f (mx 3, y 3) vào phương trình (3) sau đó giải để<br />
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)):<br />
1 x1 y1 0 0 0 1 <br />
<br />
1m 0 0 0 1 x1 y1 2 <br />
1 x y 0 0 0 <br />
me 2m 0 02 0 2 1 x y 3 Am (4)<br />
m 2 2 4<br />
3 1 x3 y3 0 0 0 5 <br />
<br />
0 0 0 1 x3 y3 6 <br />
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ử<br />
tam giác đối với trạng thái màng. Từ phương trình (4) nhận được: {} = [Am]-1 {me}<br />
Phương trình tổng quát (3) được viết như sau:<br />
1<br />
f<br />
m<br />
( x, y )<br />
m<br />
( x, y)<br />
Am me<br />
(5)<br />
3. Biểu diễn quan hệ biến dạng ( x , y ) tại điểm bất kỳ theo chuyển vị (mx , y ) và theo<br />
m<br />
<br />
<br />
<br />
chuyển vị nút me , véctơ biến dạng được viết như sau, [8]:<br />
T<br />
<br />
m<br />
( x, y ) x y xy (6)<br />
<br />
<br />
<br />
<br />
Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 27<br />
Tên lửa & Thiết bị bay<br />
<br />
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<br />
T T<br />
các biến x, y, nhận được: (mx , y ) x y xy 2 6 3 5 (7)<br />
1<br />
hoặc có thể biểu diễn như sau: (mx , y ) C m C m Am B <br />
me m me<br />
(8)<br />
ma trận biến dạng màng:<br />
y y 2 3<br />
0 y3 y1 0 y1 y2 0 <br />
1 1 <br />
m m<br />
B C A m<br />
0 x3 x 2 0 x1 x3 0 x2 x1 <br />
2 <br />
x3 x 2 y 2 y3 x1 x3 y3 y1 x2 x1 y1 y 2 <br />
m<br />
<br />
4. Quan hệ ứng suất ( x , y ) theo biến dạng và chuyển vị nút :<br />
m<br />
( x, y )<br />
me<br />
<br />
<br />
T<br />
<br />
m<br />
( x, y ) x y xy A (mx , y ) hay A Bm me<br />
m<br />
( x, y ) (9)<br />
trong đó, [A] - ma trận đàn hồi của trạng thái màng<br />
m<br />
<br />
5. Biểu diễn quan hệ ứng suất ( x , y ) với tải trọng nút tĩnh tương đương, quan hệ lực nút<br />
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à<br />
ma trận đối xứng với các thành phần kijm kij<br />
k11m k 21m k31m k41m k51m k61m <br />
m m m m m m <br />
k21 k22 k32 k42 k52 k62 <br />
1 k31 k63 <br />
m m m m m m<br />
k32 k33 k43 k53<br />
K me m m m m m m<br />
<br />
4 k41 k42 k43 k44 k54 k64 (10)<br />
k m k<br />
m<br />
k<br />
m<br />
k<br />
m<br />
k<br />
m<br />
k <br />
m<br />
<br />
51m 52<br />
m<br />
53<br />
m<br />
54<br />
m<br />
55<br />
m<br />
65<br />
m<br />
<br />
k61 k 62<br />
k 63<br />
k 64<br />
k65 k66 <br />
6. Thành lập ma trận chuyển vị ứng suất [Hme]<br />
Quan hệ ứng suất - chuyển vị: ( x , y ) H <br />
m m me<br />
(11)<br />
Trong đó: [Hm] = [Am][Bm]<br />
m<br />
Từ (11) có thể xác định được ứng suất ( x , y ) tại điểm bất kỳ (x,y) trong phần tử. Các <br />
ứ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<br />
đượ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<br />
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,<br />
y tức là ứng suất nhận được tại tâm của tam giác.<br />
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]<br />
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<br />
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<br />
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<br />
tử có 9 bậc tự do, được viết dưới dạng véctơ sau đây, [7]:<br />
T<br />
w<br />
ue<br />
1 x1 y1 w 2 x 2 y 2 w 3 x3 y3 (12)<br />
<br />
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<br />
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]:<br />
w( x, y) 1 2 x 3 y 4 x2 5 xy 6 y2 7 x3 8 ( x2 y xy2 ) 9 y3 (13)<br />
hay được viết dưới dạng ma trận sau<br />
<br />
<br />
28 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.”<br />
Nghiên cứu khoa học công nghệ<br />
<br />
w( x, y) 1 x y x2 xy y2 x3 ( x2 y xy2 ) y3 P( x, y) (14)<br />
T T<br />
trong đó, 1 2 3 4 5 6 7 8 9<br />
Đế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<br />
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.<br />
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 độ<br />
võng theo cạnh này là:<br />
wij (w) y 0 1 2 x 4 x2 7 x3 (15)<br />
w<br />
Độ dốc theo phương x hay dọc theo cạnh này là:<br />
x<br />
w w 2<br />
2 2 4 x 3 7 x (16)<br />
x ij x y 0<br />
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<br />
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<br />
w w <br />
nhiên, nếu để ý tới 4 chuyển vị nút w i , , w j , lim ta có thể thấy rằng bằng<br />
x i x j x 0<br />
<br />
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 tại hai<br />
<br />
x <br />
đ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 ;<br />
<br />
w w <br />
w tại x=0 và x=a. Từ đó, hoàn toàn có thể xác định được 4 tham số 1 ,<br />
<br />
<br />
<br />
x x <br />
i<br />
x j<br />
<br />
ij<br />
<br />
<br />
<br />
2 , 4 và 7 một cách duy nhất. Do đó, w và w hoàn toàn xác định trên cạnh biên ij và<br />
x<br />
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<br />
tử A sang phần tử B.<br />
<br />
Độ dốc theo phương y hay w dọc theo cạnh biên chung ij (y =0):<br />
y<br />
w w 2<br />
3 5 x 8 x (17)<br />
y ij y y 0<br />
<br />
Có thể thấy rằng để xác định w dọc theo cạnh biên chung ij, ta cần xác định được 3<br />
y<br />
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<br />
<br />
hai điều kiện: w w tại x = 0 và w w tại x = a<br />
<br />
y i y ij y j y ij<br />
<br />
Do đó độ dốc vuông góc với biên ij là không xác định. Độ dốc<br />
w này tại các nút i và<br />
y<br />
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<br />
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ể<br />
á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<br />
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ẽ<br />
là giá trị của w,<br />
w và w tại nút đó.<br />
x y<br />
<br />
<br />
Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 29<br />
Tên lửa & Thiết bị bay<br />
<br />
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 <br />
Đạo hàm (14) theo x và y nhận được véctơ chuyển vị, [7]:<br />
w 1 x y x 2 xy y 2 x3 ( x 2 y xy 2 ) y 3 <br />
<br />
x 0 0 1 0 x 2y 0 ( x 2 2 xy ) 3 y 2 (18)<br />
0 1 0 2 x y 0 3x 2 (2 xy y 2 ) 0 <br />
y <br />
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:<br />
ue An (19)<br />
trong đó: An - véctơ hệ số;<br />
An - ma trận hằng số<br />
1<br />
Giải phương trình (19) nhận được: Au ue<br />
(20)<br />
Đặt (20) vào (14) nhận được w = [P(x,y)] [Au]-1{ue} (21)<br />
3. Biểu diễn u<br />
( x, y ) tại điểm bất kì theo chuyển vị nút { ue<br />
}<br />
<br />
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:<br />
T<br />
T 2w 2w 2w <br />
<br />
u<br />
( x, y ) x y xy z 2<br />
y 2<br />
2 (22)<br />
x xy <br />
<br />
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:<br />
T<br />
2 2 2<br />
<br />
z x 2 P( x, y) y 2 P( x, y ) 2 xy P( x, y ) Au ue <br />
1<br />
u<br />
( x, y ) (23)<br />
<br />
Hoặc có thể viết B <br />
u<br />
( x, y )<br />
u ue<br />
<br />
<br />
Trong đó, ma trận biến dạng uốn [Bu] là [Bu] = -z[Cu][Au]-1<br />
2 <br />
2 <br />
x 0 0 0 2 0 0 6 x 2y 0<br />
2 0 0 0 0 0 2 0<br />
y 2 <br />
C u<br />
P ( x , y ) <br />
2 x 6 y (24)<br />
0 0 0 0 2 0 0 (4 x 4 y ) 0 <br />
2 <br />
2 <br />
xy <br />
4. Xác định ma trận độ cứng uốn của phần tử [Kue]<br />
Ma trận độ cứng uốn [Kue] phần tử tam giác 3 nút được xác định, [7]:<br />
1 T<br />
K ue B u <br />
Ve<br />
T t /2<br />
D Bu dV t /2 z 2 dz <br />
S<br />
A C D C A <br />
u u T u u 1<br />
dS (25)<br />
<br />
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<br />
t3<br />
C D C dS A <br />
T<br />
K ue Au <br />
12 <br />
<br />
1<br />
<br />
<br />
S<br />
u T u u 1<br />
(26)<br />
<br />
Trong đó, t và S lần lượt là chiều dày và diện tích phần tử tam giác.<br />
Công thức (26) có thể viết như sau<br />
1 T<br />
<br />
K ue Au I A u 1<br />
(27)<br />
<br />
<br />
<br />
30 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.”<br />
Nghiên cứu khoa học công nghệ<br />
<br />
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]<br />
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<br />
được ma trận có kích thước (9x9) và ma trận có tính chất đối xứng kij k ji .<br />
u u<br />
<br />
<br />
5. Thành lập ma trận chuyển vị -ứng suất [Hue]<br />
Quan hệ ứng suất - chuyển vị được cho trong phương trình (23) và (24) như sau<br />
H <br />
u<br />
( x, y)<br />
ue ue<br />
<br />
<br />
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<br />
được ứng suất tại tâm tam giác.<br />
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<br />
màng và chịu lực uốn<br />
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 độ<br />
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<br />
lực và chuyển vị tại nút 1 như sau:<br />
F1m k11m 0 0 1m <br />
<br />
<br />
F1 F1u 0<br />
<br />
k11u 0 1u <br />
<br />
(28)<br />
<br />
Tz1 0 0 0 z1 <br />
<br />
Đế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ự<br />
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<br />
diễn theo chuyển vị tại nút 3 như sau:<br />
<br />
<br />
<br />
F2m k23m <br />
0 0 3m <br />
<br />
<br />
F2 F2u 0 k23u 0 3u <br />
(29)<br />
<br />
Tz 2 0 0 0 z 3 <br />
<br />
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<br />
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<br />
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]<br />
có kích thước (18x18)<br />
k11m 0 0 k21m 0 0 k31m 0 0<br />
<br />
0 k11u 0 0 k21m 0 0 k31m 0<br />
<br />
0 0 0 0 0 0 0 0 0<br />
m m<br />
<br />
k21 0 0 k <br />
22 0 0 k32m 0 0<br />
<br />
K e 0 u<br />
k 0 0 u<br />
k 0 0 k32m 0 (30)<br />
21 22<br />
<br />
0 0 0 0 0 0 0 0 0<br />
k m 0 m<br />
0 k 0 0 k33m 0 0<br />
31 32 <br />
0 u<br />
k 0 0 u<br />
k 0 0 k33u 0 <br />
31 32<br />
<br />
0 0 0 0 0 0 0 0 0 <br />
<br />
Đố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<br />
[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<br />
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<br />
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<br />
<br />
<br />
Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 31<br />
Tên lửa & Thiết bị bay<br />
<br />
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<br />
đượ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]:<br />
- Đối với trạng thái màng<br />
AE<br />
kij kijm <br />
L<br />
- Đối với trạng thái uốn<br />
12EJ<br />
kij kiju 3 - đối với vị trí độ võng w;<br />
L<br />
6EJ<br />
kij kiju 2 - đối với vị trí góc xoay y ;<br />
L<br />
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<br />
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<br />
thân hoặc thanh dọc.<br />
2.4. Chuyển hệ trục tọa độ<br />
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ử<br />
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à<br />
đơ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<br />
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<br />
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<br />
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,<br />
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ì<br />
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ể.<br />
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<br />
dạng chuẩn đối với quan hệ lực và chuyển vị trong hệ tọa độ tổng thể như sau<br />
F *e K *e *e (31)<br />
<br />
3. MÔ PHỎNG TRƯỜNG ỨNG SUẤT - BIẾN DẠNG THÂN<br />
TÊN LỬA KH35-E SỬ DỤNG PHẦN MỀM ANSYS<br />
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<br />
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<br />
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<br />
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<br />
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<br />
đượ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ả<br />
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<br />
độ mịn của lưới chỉ đạt mức vừa<br />
phải, đáp ứng được độ hội tụ cần<br />
thiết của kết quả tính toán.<br />
Thiết lập các thuộc tính vật liệu<br />
cho các phần tử kết cấu là hợp kim<br />
nhôm AMg-6. Thiết lập kiểu liên kết<br />
giữa các bề mặt tiếp xúc dạng gắn<br />
chặt (bonded).<br />
Lựa chọn trường hợp tính toán<br />
khi tên lửa bay trong mặt phẳng<br />
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.<br />
<br />
<br />
<br />
32 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.”<br />
Nghiên cứu khoa học công nghệ<br />
<br />
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 α<br />
= 00.<br />
Tải trọng tác dụng lên thân tên lửa gồm:<br />
- Lực đẩy động cơ tác dụng lên bề mặt trong của loa phụt T = 79480 (N).<br />
- 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<br />
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).<br />
-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<br />
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<br />
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<br />
phần dọc trục theo Ox và vuông góc theo Oy tương ứng: S n x * m * g ; P n y * m * g ;<br />
Trong đó:<br />
- n x n y - hệ số quá tải theo<br />
,<br />
các phương Ox và Oy:<br />
T Fx _ kd Fy _ kd<br />
nx ; ny ;<br />
m* g m* g<br />
<br />
- m = 592,3 (kg) – khối<br />
lượng tên lửa tại thời điểm xét;<br />
- g là gia tốc trọng trường.<br />
Thiết lập tải trọng tác dụng<br />
lên ½ thân tên lửa như hình 3. Hình 3. Tải trọng tác dụng lên ½ thân tên lửa.<br />
Giải bài toán nhận được trường<br />
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.<br />
<br />
<br />
<br />
<br />
Hình 4. Trường phân bố ứng suất. Hình 5. Trường phân bố biến dạng.<br />
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<br />
= 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<br />
(σb = 3,2x108 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<br />
hạn là 2x10-3 nên kết cấu đảm bảo độ bền.<br />
4. KẾT LUẬN<br />
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<br />
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<br />
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,<br />
đã 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<br />
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<br />
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<br />
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ử.<br />
<br />
<br />
<br />
Tạp chí Nghiên cứu KH&CN quân sự, Số 40, 12 - 2015 33<br />
Tên lửa & Thiết bị bay<br />
<br />
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<br />
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.<br />
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<br />
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<br />
ở 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.<br />
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<br />
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<br />
một cách chính xác, tin cậy.<br />
TÀI LIỆU THAM KHẢO<br />
[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<br />
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).<br />
[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<br />
trong môi trường ANSYS CFX,” TC, Nghiên cứu KHCNQS, Đặc san CNTT (04-<br />
2014), tr, 92-99.<br />
[3]. Nguyễn Hoa Thịnh, Hoàng Xuân Lượng, Nguyễn Đức Cương, Trịnh Hồng Anh, Nguyễn<br />
Minh Tuấn, “Kết cấu và tính toán độ bền khí cụ bay ”, Nhà xuất bản KHKT, (2005).<br />
[4]. Лизин В.Т., Пяткин В.А, “Проектирование тонкостенных конструкций.<br />
Машиностроение”, М., (1994).<br />
[5]. Оболенский Е.П., Сахаров Б.И., Сибиряков В.А. “Прочность летательных<br />
аппаратов и их агрегатов”. Машиностроение, М., (1995).<br />
[6]. Zienkiewics O.C., “The Finite Element Methods in Engineering Science”. Mc Graw-<br />
Hill, (1979).<br />
[7]. Petyt M., “Introduction to Finite Element Vibration Analysis”. Cambridge University<br />
Press, United Kingdom, (1998).<br />
[8]. Timosenko S.P., “Theory of Plane and Shells”. Mc Graw-Hill, (1951).<br />
ABSTRACT<br />
ALGORITHM OF SETTING UP MODEL FINITE ELEMENT METHOD TO<br />
CALCULATION STRESS - DEFORMATION OF KH-35E SHIP MISLES’S HOUSING<br />
This paper presents the method and the general scheme of the model<br />
identification of aerodynamics characteristics derivatives with respect to the values<br />
of the angular velocities and accelerations of the aircraft. To exclude the confusion<br />
we call six dimensionless transient characteristics of the aircraft relatively rigidly<br />
connected with the non-deformable part of the aircraft coordinate system. The<br />
development of iterative algorithms of identification of aerodynamic characteristics<br />
and the coefficients are based on the high-precision motion model aircraft, with the<br />
possibility of monitoring the results of calculations.<br />
Keywords: Ship misles, Finite Element Methods, Theory of Plane and Shells.<br />
<br />
Nhận bài ngày 09 tháng 9 năm 2015<br />
Hoàn thiện ngày 01 tháng 12 năm 2015<br />
Chấp nhận đăng ngày 25 tháng 12 năm 2015<br />
<br />
Địa chỉ: 1Trung tâm CNCKCX - Viện KH-CNQS;<br />
2<br />
Viện KH-CNQS;<br />
*<br />
Email: binhctm@yahoo.com<br />
<br />
<br />
<br />
<br />
34 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.”<br />