Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV

NGHIÊN CỨU XÂY DỰNG HỆ PHƯƠNG TRÌNH XÁC ĐỊNH CHUYỂN VỊ
ĐÀN HỒI CỦA TAY MÁY CÔNG NGHIỆP BẰNG PHƯƠNG PHÁP PHẦN TỬ
HỮU HẠN
RESEACH ON BUILDING DEFORMATIONS ELASTIC EQUATIONS OF
INDUSTRIAL MANIPULATOR BASE ON FINITE ELEMENT METHOD
Biên Dương Xuân1, Mỳ Chu Anh1, Dũng Nguyễn Trí2
1
Học viện Kỹ thuật Quân sự, Hà Nội, Việt Nam
2
Trường Đại học Công nghiệp Hà Nội

TÓM TẮT
Khi giải quyết các bài toán cơ học và điều khiển robot, các nhà thiết kế vẫn coi các khâu
là rắn tuyệt đối và coi chuyển vị đàn hồi (CVDH) của các thành phần trong hệ thống là nhỏ,
ảnh hưởng ít đến độ chính xác hoạt động của robot. Tuy nhiên, việc nghiên cứu ảnh hưởng
của chuyển vị đàn hồi sẽ giúp cho hoạt động của tay máy chính xác hơn, thực tế hơn. Bài báo
tập trung trình bày việc xây dựng hệ phương trình chuyển vị đàn hồi và tính toán CVDH cho
hệ tay máy 2 khâu phẳng bằng phương pháp phần tử hữu hạn (PTHH). Kết quả của bài báo có
thể làm cơ sở phát triển hệ điều khiển tương ứng nhằm nâng cao độ chính xác định vị của hệ
tay máy 2 khâu phẳng nói riêng và hệ tay máy chuỗi hở nói chung khi có kể đến ảnh hưởng
của chuyển vị đàn hồi trong các khâu.
Từ khóa: robot công nghiệp, chuyển vị đàn hồi, tay máy 2 khâu phẳng.
ABSTRACT
Designers has considered links of robot are rigid bodies when they solve problems about
mechanic and control robot. They considered that deformations elastic of links are small and
don’t effect on accuracy of industrial robot. However, the study on influence of deformations
elastic will help actions of robot are more accuracy and reality. This paper focus on building
deformations elastic equations and calculation for two-link planar manipulator base on finite
element method. The results can be used to improve corresponding control system to advances

[12,13]. Ngoài ra, tác giả cũng nhắc đến phương pháp ứng dụng mô đun CAE trong các phần
mềm thiết kế CAD/CAM để tính toán chuyển vị. Tuy nhiên, phương pháp này có nhược điểm
lớn là thời gian tính toán lâu không phù hợp với điều khiển robot trong thời gian thực [14].
Một số trường hợp, cơ hệ có cấu trúc không gian nhẹ (khâu đàn hồi-flexible links) và hoạt
động với tốc độ cao thì việc xem xét ảnh hưởng của chuyển vị đàn hồi đóng một vai trò vô
cùng quan trọng trong hoạt động của cơ hệ [15,16]. Nó tác động trực tiếp tới độ chính xác
trong hoạt động.
Như vậy, vấn đề đặt ra là nếu ta xây dựng thêm được hệ phương trình xác định chuyển
vị đàn hồi của các khâu thì kết quả của nó cho phép thiết lập hệ điều khiển tương ứng chính
xác hơn, thực tế hơn do đầu vào đã kể đến ảnh hưởng của chuyển vị đàn hồi.
2. NỘI DUNG NGHIÊN CỨU
2.1. Hệ phương trình động lực học tay máy
Xét mô hình robot phẳng hai khâu [1] như Hình 1:

Hình 1. Mô hình robot phẳng hai khâu
Trong đó: Ox 0 y0 z0 , Ox1 y1 z1 , Ox 2 y2 z2 lần lượt là hệ tọa độ cố định, hệ tọa độ địa phương
gắn với khớp 1 (tại O) và khớp 2 (tại A); q1, q2 , m1 , m2 , l1 , l2 , C1 , C2 : biến khớp, khối lượng,
chiều dài khâu và vị trí trọng tâm khâu 1, khâu 2. Px , Py : ngoại lực tác dụng tại điểm thao tác
cuối B. τ 1 ,τ 2 : mô men truyền động tại các khớp 1 và 2. Gọi khoảng cách OC 1  lC 1 ,
AC 2  lC 2 ; Véc tơ biến khớp q = [ q1 , q2 ] . Ta có bảng tham số động học D-H và động lực học
T

như Bảng 1 và Bảng 2.

Bảng 1. Bảng tham số động học D-H
Khâu

Θ

d

Khâu

Vị trí trọng tâm

Ma trận mômen quán tính I

xC

yC

zC

Khối
lượng

1

l C1

0

0

m1

0

I 1y

I 1z


I yy

I zz

I xy

I yz

I zx

Ta có tọa độ điểm thao tác cuối theo hệ tọa độ cố định:
xB =l1cosq1 + l2 cos ( q1 + q2 )

(1)

yB =l1 sin q1 + l2 sin ( q1 + q2 )

Tọa độ trọng tâm các khâu theo hệ tọa độ cố định:
l1 cos q1 + lC 2 cos ( q1 + q2 ) 
lC1 cos q1 




rC1 =  lC1 sin q1  ; rC 2 = l1 sin q1 + lC 2 sin ( q1 + q2 ) 


 0 
0

τ=
( m1lC21 + m2 (l12 + lC2 2 + 2l1lC 2 cos q2 ) + I 2 z + I1z ) q1 + ( m2 (lC2 2 + l1lC 2 cos q2 ) + I 2 z ) q2
1


− 2m2l1lC 2 sin q2 q1q2 − m2l1lC 2 sin q2 q22 + (m1lC1 + m2l1 ) g cos q1 + m2lC 2 g cos(q1 + q2 )

− Px l1 sin q1 + l2 sin ( q1 + q2 )  + Py l1cosq1 + l2 cos ( q1 + q2 ) 


2
2
2
τ 2 = ( m2 (lC 2 + l1lC 2 cos q2 ) + I 2 z ) q1 + ( m2lC 2 + I 2 z ) q2 + m2l1lC 2 q1 sin q2

+ m2lC 2 g cos(q1 + q2 ) − Px l2 sin ( q1 + q2 ) + Py l2 cos ( q1 + q2 )


(4)

Dễ nhận thấy rằng, các thành phần cấu thành hệ phương trình vi phân động lực học
chưa kể đến ảnh hưởng của chuyển vị đàn hồi. Khi sử dụng hệ phương trình (4) trong giải bài
toán điều khiển sẽ xuất hiện sai số. Nếu ta xây dựng thêm hệ phương trình xác định chuyển vị
đàn hồi theo các giá trị biến khớp và có kể đến các yếu tố nội tại của các khâu như mô đun
đàn hồi (đặc trưng cho biến dạng của vật thể khi chịu tác động của ngoại lực), khối lượng
riêng (đặc trưng cho mật độ của vật chất) để đưa thêm dữ liệu cho bài toán điều khiển thì hoàn
toàn có thể làm giảm bớt sai số định vị khi tay máy hoạt động.
2.2. Xây dựng hệ phương trình xác định chuyển vị đàn hồi
Để xác định chuyển vị của các phần tử, ta có thể sử dụng một số phương pháp khác
nhau như phương pháp lực, phương pháp chuyển vị và phương pháp phần tử hữu hạn (PTHH)

{R} = {R1 , R2 , R3 , R4 , R5 , R6 , R7 , R8 , R9 }

là véc tơ chuyển vị đàn hồi tại các nút O1 , O2 , E (chuyển
vị nút) và lực nút ứng với các hệ tọa độ địa phương. (d 1 , d 4 , d 7 ) là chuyển vị dọc trục dầm
(theo trục x). (d 2 , d 5 , d 8 ) là chuyển vị theo trục y và (d 3 , d 6 , d 9 ) là chuyển vị thẳng theo trục
z. E , l1 , J1 , ρ1 , l2 , J 2 , ρ 2 lần lượt là mô đun đàn hồi, chiều dài, mô men quán tính và khối lượng
riêng

các khâu,

biến

khớp

q = [ q1 , q2 ] .
T

Gọi

{d } = {d , d , d , d , d , d , d , d , d } ,
'

'
1

'
2

'
3


'
4

'
5

'
6

'
7

'
8

'
9

hệ tọa độ cố định Ox 0 y0 z0 . Việc chuyển đổi giữa các hệ tọa độ được thực hiện qua công thức:
−1
T
=
{di } = [T ]i {di' } hay
]i {di } [T ]i {di }
{di'} [T=

(5)

Trong đó: [T ]i là ma trận chuyển đổi hệ tọa độ và trực giao [2] nên [T ]i = [T ]i

0
0
c
os(q
)
sin(
q
)
0
i
i



0
0
0 − sin(qi ) cos(q i ) 0 


0
0
0
0
0
1


cos(q i )

 sin(qi )

cos(q i ) 0 

0
1

T

(6)

Hệ phương trình cân bằng của hệ trong hệ cố định theo lý thuyết PTHH:

{R } =  K  {d }
'

'

'

(7)

Trong đó:  K '  là ma trận độ cứng của hệ trong hệ tọa độ cố định. Đưa vào các điều
kiện biên tại nút O 1 nối đất thì phương trình cân bằng có dạng mới [2]:

{R } =  K  {d }
*

*

*


 0

 0
[K ]i =  EF
−
 l
 0


 0


0

0

12 EJ
l3
6 EJ
l2

6 EJ
l2
4 EJ
l

0

0


12 EJ
l3
6 EJ
l2
0

12 EJ
l3
6 EJ
− 2
l



6 EJ 

l2 
2 EJ 
l 

0 

6 EJ 
− 2
l 
4 EJ 

l 
0


− Di Cqi2 − N i Sqi2
Ci Cqi

−Ci Sqi
Ci Cqi
Ai
Ci Sqi
−Ci Cqi
Bi

− Di Sqi2 − N i Cqi2
Di Sqi Cqi − N i Sqi Cqi
Ci Sqi
Di Sqi2 + N i Cqi2
− Di Sqi Cqi + N i Sqi Cqi
Ci Sqi

Di Sqi Cqi − N i Sqi Cqi
− Di Cqi2 − N i Sqi2
−Ci Cqi
− Di Sqi Cqi + N i Sqi Cqi
Di Cqi2 + N i Sqi2
Ci Cqi

−Ci Sqi 

Ci Cqi 
Bi 

Ci Sqi 

li

li

Mô đun đàn hồi của vật liệu khâu (coi các khâu có cùng thành phần vật liệu); J i : Mô men
quán tính của các khâu; Fi : tiết diện từng khâu. Căn cứ theo (10) ta dễ dàng xác định được
 K1'  ,  K 2'  và ta kí hiệu các phần tử của chúng như công thức (11) dưới đây:

(11)

144


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
Ma trận độ cứng toàn hệ [ K ] có dạng ma trận vuông 9x9 ứng với các chuyển vị nút.
Tuy nhiên, khi xét đến điều kiện biên ta có thể loại bỏ hàng thứ j và cột thứ j tương ứng với
chuyển vị nút q}j = 0 . Suy ra ma trận toàn hệ suy biến về dạng:

{

*
 K 44
 *
 K 54
 .
 K *  = 
 .
K*
 84
*

K85*

.
.

*
K95

*
77

K
.
.

K88*
*
K98

*

K 49
* 
K 59

. 

. 
K89* 


=
K 55'(1) + K 55'(2) ; K 66
=
K 66'(1) + K 66'(2) ; K 77
=
K 77'(2) ; K88
=
K88'(2) ; K 99
=
K 99'(2)

Tương tự, ta có thể tính được các phần tử khác của  K *  .
* Xác định véc tơ tải trọng quy nút { R*}
Lấy giả thiết vật liệu trên các khâu có tính liên tục và đẳng hướng nên có thể coi là phân
bố đều về mặt khối lượng. Gọi lực phân bố do khối lượng khâu gây ra là p (kN / cm) . Mô
hình tính tải nút có thể xét như Hình 3.

Hình 3. Mô hình tính lực nút với tải trọng phân bố
Để tính

{R } cần tính được {R} , {R} . Gọi
'

*

i

i

x là giá trị biến chiều dài có giới hạn


li

x
x2
x3
H1( x ) =
1 − ; H 2( x ) =
1− 3 2 + 2 3
l
l
l

i
pli2 i i
x 2 x3
x
=
=
τ
=
=
−
+ 2 ; H 4( x ) =
R3i =
H
p
dx
R
H

R5i =

∫H
0

5( x )

p( x ) dx =

pli i
; R6 =
2

∫H
0

6( x )

p( x ) dx = −

pli2
;
12

x2
x3
x 2 x3
H 5( x ) =
3 2 − 2 3 ; H 6( x ) =
− + 2

i
Cqi R4 − Sqi R5 
 Sq R i − Cq R i 
i 5
 i 4
i


R6



 R1'(i ) 
 '(i ) 
 R2 
 R '(i ) 
'
3 
=
{R}i =
'( i ) 
 R4 
 R '(i ) 
 5 
 R6'(i ) 

(15)

Cụ thể:


R3*

R4*

R5*

R6*}

T

(18)

Thay (12) và (18) vào (8) ta có thể tính được chuyển vị nút:

{d } =  K
*

*



−1

{R }

(19)

*

2.3. Mô phỏng kết quả tính toán chuyển vị đàn hồi

khâu B khâu
(cm) H(cm)

Chiều
dài
khâu
(cm)

Biến khớp
(rad)
0≤t ≤π

1

2,1.106

7850.10-6

0

0

3

3

20

π +t
6

=
6 2

Dễ thấy, chuyển vị tại nút O 2 lần lượt theo các trục X, Y, Z là: 0.012 mm, -0.02 mm, 0.0018 mm; tại nút E (điểm thao tác cuối): 0.012 mm, -0.04 mm, -0.0019 mm. Các giá trị này
rất nhỏ và giá trị có dấu (-) thể hiện ngược chiều so với chiều trục ứng với các nút. Vì xét tay
máy phẳng nên ta chưa xét đến chuyển vị theo chiều trục Z. Các giá trị chuyển vị này có thể
sẽ lớn nếu chịu thêm các ngoại lực tác dụng (tải trọng bên ngoài). Lúc này, vị trí điểm thao tác
cuối được mô tả như Bảng 4.
Bảng 4. Tọa độ điểm thao tác cuối
Tọa độ điểm thao tác cuối (E)

X E (mm)

YE (mm)

Khi chưa kể đến chuyển vị đàn hồi

150 3

150

150 3 + 0.012

150 − 0.04

Khi đã kể đến chuyển vị đàn hồi

* Xét vị trí tay máy với các biến khớp thay đổi theo thời gian t(s):
π
q1=

trị biến khớp thay đổi nhằm so sánh, đánh giá sai lệch vị trí điểm thao tác cuối khi có kể đến
ảnh hưởng của chuyển vị đàn hồi và làm cơ sở xây dựng các thông số đầu vào khi thiết lập hệ
điều khiển sau này.
TÀI LIỆU THAM KHẢO
[1]. Nguyễn Văn Khang, Chu Anh Mỳ, Cơ sở robot công nghiệp, Nhà xuất bản Giáo dục,
2010.
[2]. Vũ Khắc Bảy, Bài giảng phương pháp số (phương pháp phần tử hữu hạn), Bộ môn
Toán, Đại học Lâm nghiệp, 2012.
[3]. Trần Ích Thịnh, Ngô Như Khoa, Phương pháp phần tử hữu hạn, Đại học Kỹ thuật Công
nghiệp - Đại học Thái Nguyên, 2007.
[4]. Nguyễn Công Trí, Nguyễn Thị Hiền Lương, Tính toán khung phẳng bằng phương pháp
phần tử rời rạc biến thể sử dụng mô hình chuyển vị, Tạp chí KH&PT công nghệ, 2006.
[5]. Wisama Khalil, Philippe Lemoine, Autonomous calibration of robot using planar
points, International Symposium on Robotics and Manufacturing, Montpellier, France,
2009.
[6]. A. Watanabe, S. Sakakibara, A kinematic calibration method for industrial robots using
autonomous visual measurement, Annals of the CIRP vol 55/1/2006.
[7]. Ali Nahvi, John M Hollerbach, Cablration of a parallel robot using multiple kinematics
closed loops, IEEE 1994.

148


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
[8]. Andrew B Lintott, Geometric modeling and accuracy enhancement of parallel
manipulators, University of Canterbury Christchurch New Zealand, 2000.
[9]. In-chul Ha, Kinematic parameter calibration method for industrial robot manipulator
using the relative position, Journal of mechanical science and technology, 2008.
[10]. Mathieu Rognant, A systematic proceduce for the elasto-dynamic modeling and
identification of robot manipulator, IEEE Transaction on Robotics 1085-1093, 2010.