Vấn đề tiếp xúc trong việc thiết kế mô phỏng bằng phương pháp phần tử hữu hạn

VẤN ĐỀ TIẾP XÚC TRONG VIỆC THIẾT KẾ MÔ PHỎNG 
BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 
THE CONTACT PROBLEM IN THE SIMULATED DESIGN 
CREATED BY THE FINITE ELEMENT METHOD 
TRẦN QUỐC VIỆT 
Trường Cao đẳng Công nghệ, Đại học Đà Nẵng 
TÓM TẮT 
Bài báo này giới thiệu về các nghiên cứu lý thuyết xoay quanh vấn đề tiếp xúc giữa các vật 
thể, và sự áp dụng các thành tựu của những nghiên cứu này khi thiết kế mô phỏng bài toán 
có tiếp xúc bằng phần mềm ứng dụng phương pháp phần tử hữu hạn, ANSYS. 
ABSTRACT 
This article presents the theoretical research on the contact problem between two or many 
objects and the application of their results to produce the simulated design of the contact 
problems using the finite element method software ANSYS. 
1. Đặt vấn đề 
 Cùng với sự xuất hiện của máy tính và sự tiến bộ vượt bậc của kỹ thuật tính toán, 
phương pháp phần tử hữu hạn (PTHH) cũng không ngừng phát triển. Ngày nay phương pháp 
PTHH đã và đang được dùng rộng rãi trong các phần mềm mô phỏng số. 
 Trong các bài toán tính ứng suất cũng như biến dạng của các kết cấu trong cơ học chất 
rắn, nếu có từ hai vật thể trở lên tiếp xúc với nhau, thì việc tính toán cũng như việc lập mô 
hình để diễn tả các vật thể này trong các phần mềm 
áp dụng phương pháp PTHH trở nên phức tạp. Ví dụ, 
ta có một khối kim loại 1 đặt lên một khối kim loại 2, 
và khối 1 chịu tác dụng một lực F (hình 1). Khi dùng 
phương pháp PTHH để tính toán cho ví dụ này, nếu 
xét đơn giản thì có thể tính riêng khối 1: bằng cách 
chia lưới phần tử khối này với một lực tập trung F đặt 
ở mặt trên và các nút ở mặt dưới (tiếp xúc với khối 2) 
được khống chế các bậc tự do cần thiết (ví dụ chuyển 
vị theo phương thẳng đứng bằng 0). Sau đó tình toán 
khối 2: với lực phân bố do F tạo ra trong vùng tiếp 
xúc với khối 1 và các nút ở mặt dưới bị khống chế 
các bậc tự do cần thiết. Nếu mô phỏng bài toán này 
theo cách trên thì kết quả sẽ kém chính xác, bởi vì mô 
hình đó không xét đến sự tác dụng tương hỗ của các phần tử kim loại trong vùng tiếp xúc giữa 
hai khối. Thêm vào đấy, nếu khối 1 từ trên cao rơi xuống đụng vào khối 2, thì trong trường 
hợp này lại cần phải xét đến ứng xử của vật liệu khi hai khối chạm vào nhau, nếu không thì 
khối 1 sẽ đi xuyên qua khối 2 mà không “biết” đến sự hiện diện của khối 2. 
 Như vậy, việc nghiên cứu tác dụng tương hỗ của các phần tử chất rắn khi tiếp xúc với 
nhau để áp dụng vào tính toán ứng suất và biến dạng là rất cần thiết. 
2 
F 
1 
Hình 1. Ví dụ về bài toán tiếp xúc 
2. Các nghiên cứu lý thuyết về vấn đề tiếp xúc 
 Vấn đề tiếp xúc là một trong những bài toán phi tuyến phức tạp. Điểm phức tạp thứ 
nhất là do các vùng tiếp xúc luôn thay đổi không biết trước, và điểm phức tạp thứ hai là do kết 
quả tính lực ma sát trong từng giai đoạn rất là hỗn độn, do đó sẽ gây khó khăn cho việc đạt 
được một lời giải hội tụ (bằng phương pháp số). 
 Việc phân tích vấn đề tiếp xúc là cần thiết để áp dụng phương pháp PTHH trong các 
lĩnh vực xây dựng và cơ khí. Sự tiếp xúc giữa hai hoặc nhiều vật rắn biến dạng thường xảy ra 
trong các bài toán cơ học. Ngay cả trong trường hợp vật liệu là đàn hồi tuyến tính, thì vấn dề 
tiếp xúc lúc đó cũng là bài toán phi tuyến. Thực vậy, các hiện tượng tiếp xúc và ma sát được 
biểu diễn bằng các bất phương trình phi tuyến mà chúng sẽ tác động đến các chuyển vị (hoặc 
tốc độ chuyển vị) của một phần của bề mặt biên và các phản lực tiếp xúc. Các phản lực này và 
các diện tích tiếp xúc cũng là các ẩn số, chúng thay đổi dần dần khi người ta đặt một tải trọng 
ngoài. Có nhiều phương pháp đã được đề nghị để giải các bài toán tiếp xúc bằng phương pháp 
phần tử hữu hạn: phương pháp lập chương trình toán học 1, phương pháp cản trở 
(penalization) 2, phương pháp nhân thức Lagrange 3, ... 
 Feng và Touzot đã đề nghị một phương pháp hỗn hợp các phần tử hữu hạn [4]. 
Phương pháp này khá đơn giản và rất thích hợp để cho lời giải của bài toán phi tuyến cục bộ . 
Sau đây ta đi xem xét phương pháp này. 
 Nguyên tắc của phương pháp hỗn hợp: 
 Khảo sát hai vật rắn V1 và V2 tiếp xúc với 
nhau (hình 2), chịu các lực tác động P và các 
chuyển vị bị khống chế U . Khi sự tiếp xúc diễn ra, 
mỗi cặp nút tiếp xúc có thể rơi vào ba trường hợp: 
tiếp xúc dính, tiếp xúc trượt và không tiếp xúc. 
Trong ba trường hợp này, các điều kiện tiếp xúc và 
các phương trình cân bằng các lực tiếp xúc được 
trình bày trong [4]. Trong bối cảnh của các phần tử 
hữu hạn, phương trình cân bằng của hai vật rắn 
được viết trực tiếp như sau: 
      K U P R  (1) 
trong đó: 
K: ma trận độ cứng của cấu trúc 
U: véc tơ chuyển vị của cấu trúc, P: véc tơ 
ngoại lực, R: véc tơ phản lực tiếp xúc. 
 Vì U và R chưa biết, phương trình (1) không thể được giải trực tiếp. Véc tơ R 
trước tiên được xác định cục bộ bằng phương pháp lực. Tiếp theo, véc tơ U được tính toán 
toàn bộ bằng phương pháp chuyển vị. Vì vậy, phương pháp này được gọi là “phương pháp 
hỗn hợp”. 
 Phương trình (1) được viết lại dưới dạng sau: 
* *
0rr rc rr
T
rc cc c
K K PU
K K PU R
      
       
      
 (2) 
trong đó U* tương ứng với các nút có tiếp xúc, Ur liên quan đến các nút còn lại trong kết 
cấu, R* là véc tơ phản lực liên quan đến các nút tiếp xúc. Bằng cách khử Ur, phương 
trình (2) trở thành: 
U
P 
P 
V1 
V2 
R 
R 
Hình 2. Kết cấu đang tiếp xúc 
      * * * *cc cK U P R     (3) 
trong đó        
1* T
cc cc rc rr rcK K K K K

     (4) 
          1* Tc c rc rr rP P K K P

  (5) 
Theo phương trình (3), ta có được: 
      * * * *c pU F R U    (6) 
với 
1
* *
c ccF K

       : ma trận mềm dẻo (matrix of flexibility), thu hẹp cho các nút tiếp xúc. 
    * * *p c cU F P    : véc tơ chuyển vị gây ra bởi ngoại lực. 
- Phương trình độ cách giữa các cặp nút tiếp xúc: 
Người ta lập độ cách giữa các cặp nút tiếp xúc khởi đầu từ các chuyển vị của các nút 
tiếp xúc và các độ cách ban đầu  *0X giữa hai vật rắn V1 và V2: 
        * * * *1 2 0X U U X   (7) 
Bằng cách phối hợp các phương trình (2) và (4), ta có: 
        * * * * *0pX F R U X     (8) 
trong đó: 
  *pU : véc tơ các độ cách gây ra bởi ngoại lực, 
 * * *1 2c cF F F            : ma trận mềm dẻo. 
- Quan hệ trong hệ thống tọa độ cục bộ: 
 Để thỏa mãn các điều kiện tiếp xúc và thiết lập các phương trình cân bằng của tiếp 
xúc, người ta chuyển phương trình (8) vào hệ tọa độ cục bộ n – t, với n là phương pháp tuyến 
và t là phương tiếp tuyến với bề mặt tiếp xúc. Giả thiết rằng x và r lần lượt là véc tơ các 
độ cách và véc tơ các lực tiếp xúc trong hệ tọa độ cục bộ thì ta có: 
 x= TX* (9) 
 r= TTR* (10) 
trong đó T là ma trận xoay tương ứng với phép biến đổi tọa độ giữa hệ tọa độ cục bộ và hệ 
tọa độ tổng thể. 
Thay biểu thức (9) và (10) vào phương trình (8), ta được phương trình độ cách giữa 
các cặp nút tiếp xúc trong hệ tọa độ cục bộ: 
 x= fr + up + x0 (11) 
trong đó 
 f = TF*TT, up= TUp
*, x0= TX0
* 
- Dạng số gia (increment form) của phương trình độ cách: 
Trong các bài toán tiếp xúc, ma sát sinh ra do tác dụng tương hỗ giữa hai vật thể tại 
các bề mặt tiếp xúc của chúng. Dưới tác dụng của tải, hai vật thể có thể có chuyển động tương 
đối với nhau. Ta biết rằng chuyển động này phụ thuộc vào các phản lực pháp tuyến và tiếp 
tuyến trên bề mặt tiếp xúc. Nó xuất hiện ở một ngưỡng nào đó. Khi vượt qua ngưỡng này sẽ 
xuất hiện sự trượt. Ngoài ra các hiêu ứng ma sát còn phụ thuộc vào sự tiến triển theo thời gian 
của tải. Do vậy, một phép xấp xỉ theo kiểu số gia là cần thiết để theo sát hơn sự tiến hóa của 
sự trượt. Người ta đã áp dụng quá trình từng bước như sau: 
Dùng ký hiệu “” để diễn tả các biến số gia. Ở bước (i), phương trình (11) thành: 
 xi= fri + upi + x0 (12) 
trong đó 
 ri= ri-1 + ri, upi= upi-1 + upi (13) 
 Thay (13) vào (12), ta có: 
 xi= fri + upi + fri-1 + upi-1 + x0 (14) 
hơn nữa: xi-1= fri-1 + upi-1 + x0 (15) 
 Từ đó ta sẽ nhận được dạng số gia của phương trình (11) như sau: 
 xi = f ri + upi + xi-1 (16) 
 Khi kết thúc bước (i-1), xi-1 đã biết. Ở bước i, số gia của tải trọng ngoài Pi đã 
cho, upi có thể được tính toán bằng cách giải hệ thống tuyến tính Kupi = Pi. 
- Phương trình dính: 
Ban đầu, ta giả thiết rằng tất cả các nút tiếp xúc đều ở trạng thái dính, hay nói khác đi 
là độ cách bằng 0: 
 xi = 0 (17) 
 Vậy phương trình dính sẽ là: 
 fri = b (18) 
trong đó 
 b = - upi - xi-1 
- Sửa đổi ma trận mềm dẻo và thành phần thứ hai (b): 
Trong trường hợp tiếp xúc dính, người ta giải trực tiếp phương trình (18) không có sự 
sửa đổi nào. Phương trình này cũng luôn còn giá trị cho các nút trong trạng thái trượt hoặc 
không tiếp xúc, nhưng trong hai trường hợp này, f và b phải được thay đổi (4 ). 
- Tiếp xúc 3D: 
Các phương trình căn bản của bài toán tiếp xúc 3D cũng giống như trên. Tuy nhiên 
trong trường hợp 3D, sự thay đổi ma sát sẽ gây khó khăn hơn. Trong bài toán 2D, sự trượt 
hoặc sang trái hoặc sang phải, do đó việc xác định hướng trượt D tương đối dễ. Nhưng trong 
bài toán 3D, việc tìm ra hướng trượt rất khó khăn. Để giải quyết vấn đề này ta có thể dùng 
phương pháp chiếu (projection) 5, hoặc phương pháp tường minh (explicit method) 4. 
3. Áp dụng lý thuyết tiếp xúc để thiết kế mô phỏng trong phần mềm ANSYS 
Khi thiết kế mô phỏng một bài toán 
cơ học có tiếp xúc (ví dụ vật 1 đi xuống 
tiếp xúc với vật 2, hình 3) trong phần mềm 
ANSYS, thì ngoài việc chọn kiểu phần tử 
cho vật 1 và vật 2, ta còn phải chọn kiểu 
phần tử tiếp xúc cho các nút có khả năng 
tiếp xúc khi hai vật chạm nhau. Nhìn vào 
hình 3, ta có thể tiên đoán ngay rằng các 
nút có khả năng tiếp xúc là các nút ở cạnh 
Vật 1, kiểu phần tử 1 
 Vật 2, kiểu 
 phần tử 2 
Hình 3. Bài toán mô phỏng 2 vật tiếp xúc 
trái và cạnh dưới của vật 1 và các nút ở cạnh trên của vật 2. Như vậy, ta phải chọn một kiểu 
phần tử tiếp xúc nào đó đã được định nghĩa trong ANSYS để gán cho các nút trên, sao cho 
ứng xử của các nút này khi tiếp xúc phù hợp với sự tiếp xúc trong thực tế nhất. 
Áp dụng triệt để các thành tựu của các nghiên cứu lý thuyết về tiếp xúc, ANSYS đã 
đưa ra 3 mô hình tiếp xúc: nút - nút, nút - bề mặt, và bề mặt - bề mặt 6. Trên cơ sở đó ta có 
thể gặp các kiểu phần tử tiếp xúc trong ANSYS như sau: 
- Các kiểu phần tử tiếp xúc “nút - nút”: các kiểu phần tử CONTAC12 và 
CONTAC52 sử dụng khi có sự trượt nhỏ và ma sát theo định luật Coulomb. 
- Các kiểu phần tử tiếp xúc “nút - bề mặt": các kiểu phần tử CONTAC26 và 
CONTAC48 sử dụng cho bài toán 2D và CONTAC49 cho 3D. Các kiểu CONTAC48 và 
CONTAC49 ứng dụng khi có sự trượt lớn, sự biến dạng lớn và sự lưới hoá khác nhau giữa 
các thành phần tiếp xúc. Kiểu CONTAC26 dùng để lập mô hình tiếp xúc nút mềm với bề mặt 
cứng. 
- Các kiểu phần tử tiếp xúc “bề mặt - bề mặt”: ANSYS đưa ra các phần tử tiếp xúc 
“bề mặt cứng với bề mặt mềm”. Bề mặt cứng được gọi là “mục tiêu” (target) và được mô hình 
hoá bằng kiểu phần tử TARRGE169 cho trường hợp 2D và TARGE170 cho 3D. Bề mặt của 
vật có thể biến dạng (bề mặt mềm) gọi là mặt “tiếp xúc” (contact) được mô hình hóa bằng 
các kiểu phần tử CONTA171, CONTA172, CONTA173 hoặc CONTA174. 
Tổng quát là như vậy, nhưng để có được một kết quả hội tụ khi thiết kế mô phỏng 
trong ANSYS, thì người sử dụng phải nghiên cứu kỹ các kiểu phần tử trên, trong bài báo này 
chúng tôi chưa thể trình bày được. 
4. Kết luận 
 Lý thuyết tiếp xúc đã được nhiều nhà khoa học quan tâm nghiên cứu. Những thành tựu 
về lý thuyết này cũng đã được sử dụng trong các phần mềm áp dụng phần tử hữu hạn để thiết 
kế mô phỏng các bài toán cơ học có tiếp xúc, trong đó một phần mềm rất mạnh là ANSYS đã 
áp dụng rất triệt để các thành tựu lý thuyết này. 
TÀI LIỆU THAM KHẢO 
[1] Feng Z. Q., Résolution du problème de contact unilatéral par une méthode de 
programmation mathématique: LCP – Linear Complementarity Problem, rapport 
interne, MNM/UTC, 1990. 
[2] De Saxcé G., Feng Z. Q., Touzot G., Rigid – plastic implicit schema for tow and three 
dimensional analysis of metal forming by finite element method, Engrg. Comp., 1992. 
[3] Alart P., Curnier. A., A mixed formulation for frictional contact problems prone to 
Newton like methods, Comp. Meth. Appli. Mech. Engng., 1992. 
[4] Feng Z. Q., Touzot G., Analyses bi – et tridimensionnelle de problèmes de contact 
avec frottement par une méthode mixte des éléments finis, Revue eueurropes enne des 
éléments finis, Vol.1 - n
o
 4, pp 441 – 459, 1992. 
[5] Bendhia H., Durville D., Two-dimensional modeling of contact-friction phenomena in 
the blankholder arearrea for the drawing process, Euromech 273, Unilateral contact 
and dry friction, Montpellier, France, 1990. 
[6] ANSYS, Structural Analysis Guide, 3rd Edition., SAS IP Inc., USA, 1998.