Luận án Nghiên cứu tính toán tác dụng xuyên của bom đạn vào kết cấu công sự bằng phương pháp phần tử rời rạc

pháp tuyến và tiếp tuyến; 
 61 
 Dn: lực cản tiếp xúc theo phương pháp tuyến; 
 Ds: lực cản tiếp xúc theo phương tiếp tuyến; 
  1n n n xi xj yi yj
N
2
D c n c V V cos V V sin
  (2.64) 
xi xj yi yj
s s s
1i i j j
N
2
V V sin V V cos
D c s c
d d
  
 
   

 
 (2.65) 
 Dni = - Dnj = - Dn và Dsi = - Dsj = Ds (2.66) 
 Trong hệ tọa độ tổng quát: 
 xi ni siD D cos D sin và yi ni siD D sin D cos (2.67) 
 xj nj sjD D cos D sin và yj nj sjD D sin D cos (2.68) 
 Hệ số cản luôn luôn tỉ lệ thuận với hệ số cản giới hạn: 
 critn n nc c  và 
crit
s s sc c  (2.69) 
Hệ số cản giới hạn: criti i
i
m
c 2m 2
k
  (i=n,s) (2.70) 
trong đó: 
βn,s : hệ số tỉ lệ được xác định bằng thực nghiệm; 
 i : tần số dao động riêng của hệ khi không có cản; 
 ki : độ cứng lò xo tiếp tuyến của phần tử; 
 m : khối lượng hiệu quả của hệ, 1 2
1 2
m m
m
m m
; 
 m1 và m2 lần lượt là khối lượng của phần tử 1 và 2. Trong trường hợp 
tiếp xúc giữa phần tử và biên thì m là khối lượng phần tử. 
2.5.3.2. Cản toàn bộ 
 Lực cản toàn bộ (global damp) là tập hợp các lò xo giảm chấn giữa phần 
tử và môi trường xung quanh. Lực cản toàn bộ tỉ lệ với giá trị tuyệt đối của vận 
tốc chuyển động của phần tử và có chiều ngược lại với véc tơ vận tốc.
 62 
    dF C U  và    cC M (2.71) 
 c : hệ số cản được xác định bằng thực nghiệm. 
 Các mô hình trong PFC [47]: d d iF F sign V , trong đó: (2.72) 
 d: hệ số cản được xác định bằng thực nghiệm, tương tự c; 
 Fi: thành phần lực tổng quát, bao gồm: lực tiếp xúc, lực bản thân và 
ngoại lực tác dụng lên phần tử; 
 V: vận tốc chuyển động của phần tử; 
1
sign(V) 1
0
 nếu 
V 0
V 0
V 0
 
 
 
 (2.73) 
2.5.4. Mô hình liên kết (bond model) 
 Các phần tử có thể liên kết với nhau tại điểm tiếp xúc, có hai loại liên 
kết: liên kết đơn (contact bond) và liên kết song song (parallel bond). Liên kết 
đơn, liên kết hai phần tử xung quanh điểm tiếp xúc, trong khi liên kết song 
song liên kết hai phần tử ở hai bên điểm tiếp xúc. Phần tử liên kết có dạng 
thanh phẳng, có tiết diện là hình chữ nhật nếu hai phần tử tiếp xúc là hình đĩa 
tròn (disk), tiết diện hình tròn nếu hai phần tử tiếp xúc là hình cầu (sphere). 
Phần tử liên kết đơn tồn tại 2 thành phần lực (Fx,Fy), trong khi phần tử liên kết 
song song tồn tại 2 thành phần lực và mô men (Fx,Fy,M). Trong DEM chỉ có 
liên kết giữa hai phần tử, không tồn tại liên kết giữa phần tử và tường. 
Fy
Fx
Fy
M
X0
Y
(a) (b)
Fx
Hình 2.17 Liên kết đơn (a) và liên kết song song (b) 
 63 
2.5.4.1. Liên kết đơn 
 Liên kết đơn giữa hai phần tử giống như một cặp lò xo đàn hồi theo 
phương pháp tuyến và tiếp tuyến. Lực pháp tuyến và tiếp tuyến của lò xo bị 
giới hạn bởi cường độ vật liệu chịu kéo và cắt của phần tử liên kết. Các thành 
phần lực đó gia tăng khi hai phần tử có xu hướng tách nhau ra (Un<0), trong 
trường hợp này phần tử liên kết có tác dụng kéo hai phần tử trở lại. Nếu như 
lực kéo tiếp xúc pháp tuyến đạt đến giới hạn cường độ chịu kéo của vật liệu thì 
liên kết sẽ bị phá vỡ, lực tiếp xúc pháp tuyến và tiếp tuyến sẽ đặt lại bằng 
không. Nếu như lực cắt tiếp xúc đạt đến giới hạn cường độ chịu cắt của vật liệu 
thì liên kết sẽ bị phá vỡ, lực tiếp xúc không thay đổi và trở thành lực nén. 
Un
Fn Fs
Us
Fnc
tr­ît xuyªn
Kn
1
ks
Fsc
Fsmax
ph¸ vì liªn kÕt
tr­ît
ph¸ vì liªn kÕt
(khi Un>0)
(a) (b)
1
Hình 2.18 Thành phần lực tiếp xúc pháp tuyến (a) và tiếp tuyến (b) 
 Mối quan hệ giữa lực và chuyển của phần tử liên quan đến các thành 
phần lực tiếp xúc pháp tuyến và tiếp tuyến tại điểm tiếp xúc, các liên kết tiếp 
xúc và ứng xử trượt xảy ra đồng thời. Trong hình vẽ trên, Fn là lực pháp tuyến 
(khi Fn>0 là lực kéo), Un là chuyển vị pháp tuyến (khi Un>0 hai phần tử có tiếp 
xúc với nhau), Fs là độ lớn của tổng lực tiếp xúc tiếp tuyến và Us là độ lớn của 
tổng chuyển vị tiếp tuyến. Phần tử liên kết là thanh phẳng có tiết diện chữ nhật 
và chiều dài nối tâm hai phần tử, có các tham số sau: 
 64 
 - Bán kính : i j
1
R r r
2
  , m; (2.74) 
 - Chiều dài : i jL 2R r r 
 , m; (2.75) 
 - Diện tích mặt cắt ngang : A 2RB  , m2; (2.76) 
 - Mô men quán tính góc : 3
2
I BR
3
  , m4; (2.77) 
 - Độ cứng tiếp xúc pháp tuyến : cn
AE
k
L
 , N/m; (2.78) 
 - Độ cứng tiếp xúc tiếp tuyến : cs 3
IE
k 12
L
 , N/m; (2.79) 
 Ở cấp độ vi mô (micro level), hệ số nở ngang  không ảnh hưởng đến độ 
cứng pháp tuyến và tiếp tuyến. Với các mô hình cao hơn (macro level) thì hệ 
số  ảnh hưởng trực tiếp đến tỉ lệ n sk / k giữa độ cứng pháp tuyến và tiếp 
tuyến. Thực tế, sau khi đã xác định được các tham số Ec và kn, hệ số ks được 
xác định bằng công thức sau khi hiệu chỉnh: 
n
s
n s
k
k
k / k
 (2.80) 
 Lực kéo pháp tuyến và tiếp tuyến lớn nhất: 
 nc cF 2RB   và sc cF 2RB   (2.81) 
 Như vậy, ở cấp độ vi mô liên kết đơn bao gồm 5 tham số: 
 Ec, n sk / k , µ, c và c (2.82) 
2.5.4.2. Liên kết song song 
Liên kết song song giữa hai phần tử được đặc trưng bằng cặp lò xo đàn 
hồi theo phương pháp tuyến và tiếp tuyến. Độ lớn của lực liên kết gia tăng khi 
hai phần tử bắt đầu tách nhau ra, liên kết có thể bị phá vỡ theo phương pháp 
tuyến nếu như ứng suất chịu kéo vượt quá cường độ chịu kéo của vật liệu hoặc 
bị phá vỡ theo phương tiếp tuyến nếu ứng suất cắt vượt quá cường độ chịu cắt 
của vật liệu liên kết. 
 65 
(a) (b)
C
2R
Fni
Fsi
M3
Fsi
M3
Fni 2R
B
L
i j
L
Hình 2.19 Liên kết song song có mặt cắt hình tròn (a) và chữ nhật (b) 
Các tham số của liên kết như sau: 
- Độ cứng pháp tuyến : 
c
n
E
k
L
 , N/m3; 
- Độ cứng tiếp tuyến : 
c
s 3
IE
k 12
AL
 , N/m3; 
- Cường độ chịu kéo : c , N/m
2; 
- Cường độ chịu cắt : c , N/m
2; 
- Bán kính liên kết : i jR min r , r  , m với: 0 1  . 
trong đó:  là hệ số bán kính liên kết, thể hiện phạm vi liên kết giữa hai 
phần tử, r là bán kính nhỏ hơn của 2 phần tử. 
Liên kết song song 
 66 
Ở cấp độ vi mô (micro level), hệ số nở ngang  không ảnh hưởng đến 
các tham số phân tích. Ở cấp độ vĩ mô (macro level)  ảnh hưởng rất nhiều đến 
tỉ lệ n sk / k . Chính vì vậy, sau khi đã xác định được các tham số: cE và nk , 
hệ số sk được xác định sau khi hiệu chỉnh vật liệu: 
n
s
n s
k
k
k / k
 (2.83) 
Liên kết song song bao gồm các thành phần lực: nF , sF và 3M . 
 Véc tơ lực tổng quát: i in isF F F 
 (2.84) 
 Khi liên kết được thiết lập thì đặt: iF 0 và 3M 0 . 
 Số gia về lực và mô men trong mỗi bước thời gian: 
 ni n n iF k A U n và si s siF k A U , với: i iU V t (2.85) 
 3 n 3M k I  với: 3 3 j 3i t    (2.86)
 trong đó: 
A (m2) là diện tích mặt cắt ngang và I (m4) là mô men quán tính của 
phần tử liên kết. 
 A 2RB và 
32
I R B
3
 : mặt cắt dạng hình chữ nhật; (2.87) 
2
A R và 
41
I R
4
 : mặt cắt dạng hình tròn; (2.88) 
 Lực pháp tuyến và tiếp tuyến được cập nhật sau mỗi bước thời gian: 
 ni ni ni
N N 1
F F F
 và si si si
N N 1
F F F
 (2.89) 
 3 3 3
N N 1
M M M
 (2.90) 
 Ứng suất kéo và cắt lớn nhất được xác định từ lý thuyết về dầm: 
3ni
max
MF
R
A I
 và 
si
max
F
A
 (2.91)
 67 
 Nếu như ứng suất kéo (cắt) vượt quá cường độ chịu kéo (cắt) của vật 
liệu: max c  hoặc max c  thì liên kết sẽ bị phá vỡ. Như vậy, ở cấp độ vi 
mô liên kết song song bao gồm 5 tham số: cE , n sk / k , µ, c và c (2.92) 
2.5.5. Hiệu chỉnh vật liệu 
Khác với mô hình liên tục, các tham số phân tích được lấy theo giá trị 
thực tế khi thí nghiệm, bao gồm: mô đun đàn hồi (E), hệ số ma sát (µ), cường 
độ chịu nén của vật liệu (), cường độ chịu kéo của vật liệu (). Các tham số 
phân tích trong DEM là các tham số vi mô, bao gồm [47], [61], [62]: 
- Đối với phần tử : Ec, kn, n sk / k và . 
 - Đối với liên kết tiếp xúc : Ec, n sk / k , c và c. 
 - Đối với liên kết song song : c n n s c c, E , k , k / k , ,   . 
Như vậy, các tham số phân tích là các hàm số phụ thuộc vào giá trị vật 
liệu (E và ), do đó phải thông qua hiệu chỉnh để kết quả đạt được gần giống 
theo mô hình vật liệu thật, đó là quan hệ giữa ứng suất và biến dạng: -. Để 
làm được việc này, thực hiện thí nghiệm có nở hông: sử dụng tường (wall) di 
chuyển với vận tốc nhỏ (V), thay đổi dần dần các tham số phân tích của DEM 
(Hình b và c) cho đến khi kết quả đường cong - gần đúng với sơ đồ nén kết 
cấu thực (Hình a). 
kÕt cÊu
t­êng
t­êng
t­êng
t­êng
(a) (b) (c)
V V
V
V
phÇn tö
liªn kÕt
t­êng
V
t­êng
Hình 2.20 Mô hình kết cấu (a), sơ đồ chịu nén (b) và sơ đồ chịu cắt (c) 
 68 
Biến dạng theo các phương: 0
0
L L
2
L L
 
 (2.93) 
Vận tốc của tường: w ycu G G    (2.94) 
Số gia lực tiếp xúc của tường gia tăng khi tường chuyển động theo từng 
bước thời gian: w nw c wF k N u t  (2.95) 
Số gia ứng suất trong tường: nw c ww
k N u
t
A
  

 (2.96) 
trong đó: 
L và L0: chiều dài hiện tại và ban đầu của mẫu vật liệu; 
Nc: tổng số tiếp xúc với tường; 
knw: là độ cứng trung bình của tất cả các tiếp xúc; 
A: diện tích của tường; 
 và yc: ứng suất đo được và ứng suất giới hạn yêu cầu; 
Để hội tụ, giá trị tuyệt đối của số gia ứng suất trong tường phải nhỏ hơn 
giá trị cho phép: 
w  
 
, trong đó: là hệ số thực nghiệm, =0,05 (2.97) 
nw ck N G t
A
 
  và 
nw c
A
G
k N t
 (2.98) 
Trong phân tích DEM: giá trị độ cứng kn và ks của tường được lấy cao 
hơn nhiều so với phần tử để khẳng định không tồn tại liên kết giữa tường với 
các phần tử có tiếp xúc. Thực hiện hiệu chỉnh theo từng bước, kết quả từ thí 
nghiệm trên là xác định được tỉ số giữa độ cứng pháp tuyến và tiếp tuyến: 
( n sk / k ) và ( n sk / k ). Như vậy, trước hết tính toán trước các giá trị độ cứng 
pháp tuyến, giá trị độ cứng tiếp tuyến được xác định từ mô hình hiệu chỉnh: 
n
s
n s
k
k
k / k
 và 
n
s
n s
k
k
k / k
 (2.99) 
 69 
Để xác định mô đun đàn hồi (E) và hệ số nở ngang () của sử dụng 
thuật toán tăng tải và giảm tải, cần xác định các tham số sau: 
a: ứng suất trục; 
d: ứng suất lệch trục; 
a: biến dạng dọc trục; 
v: biến dạng thể tích; 
v=x+y (2.100) 
 Khi số gia ứng suất trục bằng ứng suất lệnh trục d a   thì: 
 a d
a a
E
  
  
 (So sánh với E thực của vật liệu) (2.101) 
vx
a a
1
  
 
  
 (So sánh với  thực của vật liệu) (2.102) 
2.5.6. Lựa chọn bán kính phần tử 
Bán kính phần tử (r) trong DEM có mối liên hệ chặt chẽ với hệ số độ 
cứng phá hủy KIc trong Cơ học phá hủy. Thực tế có 3 mode phá hủy: mode I, 
mode II và mode III. Đối với mode I [9], [15], [23], [35], [56], [63]: 
cIcK r  (2.103) 
trong đó: 
: hệ số kể đến ảnh hưởng làm giảm yếu mô men uốn (<1); 
α: hệ số kể đến ảnh hưởng do không đồng nhất của vật liệu (α 1); 
c : cường độ chịu kéo của liên kết song song; 
Như vậy trong phân tích DEM, bán kính phần tử được chọn sao cho thỏa 
mãn với cơ chế phá hủy thực của vật liệu. 
2.6. Phương trình chuyển động tổng quát của hệ 
2.6.1. Thiết lập phương trình chuyển động 
       M U C U F   , trong đó: (2.104) 
 70 
 M : ma trận khối lượng của hệ; 
  C : ma trận cản toàn bộ của hệ; 
 U : véc tơ gia tốc của hệ; 
 U : véc tơ vận tốc của hệ; 
 U : véc tơ chuyển vị của hệ; 
 F : véc tơ lực tổng hợp của hệ, bao gồm: lực tiếp xúc, ngoại lực tác 
dụng, trọng lực bản thân và lực cản cục bộ. Biểu diễn các đại lượng: 
  
T
x1 y1 1 x2 y2 2 xn yn nU U U U U ... U U    (2.105) 
  
T
x1 y1 1 x2 y2 2 xn yn nU U U U U ... U U   
         (2.106) 
  
T
x1 y1 1 x2 y2 2 xn yn nU U U U U ... U U   
         (2.107) 
  
T
x1 y1 1 x2 y2 2 xn yn nF F F M F F M ... F F M (2.108) 
  
1
1
1
2
2
2
n
n
n
m
m
I 0
m
m
M I
...
...
0 m
m
I
 (2.109) 
 F là thành phần véc tơ lực thay đổi theo từng bước thời gian, phụ thuộc 
vào trạng thái tiếp xúc của các phần tử trong hệ. Do đó, các tham số cần xác 
định cho mỗi phần tử trong từng bước thời gian bao gồm: 
 71 
 - Véc tơ chuyển vị :  
T
i xi yi iU U U  , i=1÷n (2.110) 
 - Véc tơ vận tốc :  
T
i xi yi iU U U 
   , i=1÷n (2.111) 
 - Véc tơ gia tốc :  
T
i xi yi iU U U 
   , i=1÷n (2.112) 
2.6.2. Phương pháp số giải phương trình chuyển động 
 Để giải phương trình chuyển động tổng quát của hệ, sử dụng phương 
pháp sai phân trung tâm, CDM (Central Difference Method) [46]. Sai phân 
trung tâm là giá trị trung bình của sai phân tiến và sai phân lùi. Sai phân trung 
tâm sử dụng giá trị của hai điểm đối xứng so với điểm i. Chia miền thời gian t 
đã xác định thành các phân đoạn, tại mỗi đầu phân đoạn ta có các nút sai phân. 
Tập hợp các nút sai phân tạo thành lưới sai phân, viết phương trình sai phân 
cho mỗi nút lưới ta có n phương trình tương ứng với n ẩn số (sau khi đã áp 
dụng điều kiện biên), giải hệ phương trình đại số tuyến tính này ta tìm được giá 
trị của chuyển vị tại mỗi một nút lưới. 
ti-1 ti ti+1 t
Ui-1
Ui+1
Ui
 t 
Hình 2.21 Sơ đồ sai phân trung tâm 
 Vận tốc và gia tốc được xác định ở nửa bước thời gian: 
 
  i 1 i
1
i
2
U U
U
t
 và 
  i i 1
1
i
2
U U
U
t
 (2.113) 
 
  1 1
i i
2 2
i
U U
U
t
 
 (2.114) 
 72 
 Thay thế phương trình (2.113) vào phương trình (2.114) ta được: 
 
   i 1 i i 1
2i
U 2 U U
U
t
 (2.115) 
 Để tính tích phân ta chia thành n khoảng thời gian, mỗi khoảng thời gian 
là t . Giả thiết trong mỗi khoảng thời gian từ i 1t đến i 1t các đại lượng ma 
trận khối lượng M và cản nhớt C không thay đổi, còn véc tơ tải trọng của phần 
tử thay đổi theo quy luật cho trước và không phụ thuộc vào chuyển vị U. Như 
vậy ta đã tuyến tính hóa bài toán trong khoảng thời gian t . Khi đó ta có thể 
áp dụng sơ đồ sai phân trung tâm cho bài toán tuyến tính. Giả thiết: 
Tại thời gian i 1t :  i 1U U ;  i 1U U   và  i 1U U   (2.116) 
 Tại thời gian it :  iU U ;  iU U   và  iU U   (2.117) 
 Tại thời gian i 1t :  i 1U U ;  i 1U U   và  i 1U U   (2.118) 
 Theo sơ đồ sai phân trung tâm: 
    2 i 1 i i 1i
1
U U 2 U U
t 
 (2.119) 
   i 1 i 1i
1
U U U
2 t 
 (2.120)
 Thay vào phương trình (2.104) tại thời điểm i ta có: 
       ii iM U C U F(t)   (2.121) 
  
   
 
  
 i 1 i i 1 i 1 i 12 i
U 2 U U U U
M C F(t)
2 tt
 (2.122) 
 Khai triển phương trình (2.122) ta được: 
  i 1 iM U F(t) với    2
1 1
M M C
2 tt
 (2.123)
          2 2i i i 1i
2 1 1
F(t) F(t) M U M C U
2 tt t 
 (2.124) 
 73 
 Vấn đề đặt ra là cần tìm i 1U . Điều này có thể thực hiện được bằng 
cách dựa vào điều kiện biên ban đầu: 
 Khi t=0 thì  0U U và  0U U   (2.125) 
 Thay vào phương trình (2.104) ta có: 
 
    
 
0 0
0
F C U
U
M

 (2.126) 
 Thay i 1U từ phương trình (2.120) vào phương trình (2.119) ta được: 
 
    i 1 i i 1i
2i
U 2 U 2 t U U
U
t

 (2.127) 
 
   i 1 i i
2i
U U t U
U 2
t

 (2.128) 
    2i 1 i i i
1
U U t U t U
2 
   (2.129) 
 Tại thời điểm ban đầu khi t=0 thì  1U được xác định theo công thức: 
    21 0 0 0U U t U 0,5 t U   (2.130) 
2.6.3. Xác định bước thời gian tích phân 
 Phương pháp tích phân số sẽ ổn định nếu chọn bước thời gian t nhỏ 
hơn bước thời gian giới hạn critt : 
 crit
i 1 n i 1 n
2 m
t min min 2
k  
 
 (2.131) 
 Nếu kể đến tương tác của nhiều phần tử thì: 
 crit t
i 1 n
m
t min 2
k
 
, trong đó: t 0,1÷0,15 (2.132) 
Tính toán cho thấy, khi kể đến hệ số αt thì phương pháp tích phân số 
luôn cho kết quả ổn định và chính xác. 
 74 
2.7. Thuật toán 
2.7.1. Hiệu chỉnh vật liệu 
NhËp c¸c tham sè vËt liÖu
Ch­¬ng tr×nh PFC2D
1. S¬ ®å kÕt cÊu
2. Ph¸t sinh phÇn tö
3. M« h×nh liªn kÕt
4. Tham sè cÊp vi m«
5. Biªn
6. Ch¹y ch­¬ng tr×nh
7. KÕt cÊu c©n b»ng
M« pháng PFC2D
So s¸nh kÕt qu¶ gi÷a
m« pháng vµ thÝ nghiÖmKÕt qu¶ thÝ nghiÖm
Thùc hiÖn ch­¬ng tr×nh
§óng
Sai
HiÖu chØnh vËt liÖu
Hình 2.22 Hiệu chỉnh vật liệu 
Bao gồm 5 bước: 
Bước 1: Phát sinh phần tử, các phần tử có bán kính: [rmin rmax]. 
Bước 2: Xác định độ lệch ứng suất cho phép của tường để hệ cân bằng, 
thông thường nhỏ hơn 1% ứng suất nén nở hông của mô hình. 
Bước 3: Xóa đi các phần tử có ít hơn 3 tiếp xúc. 
 75 
Bước 4: Gán liên kết song song giữa các phần tử có tiếp xúc với nhau. 
Bước 5: Xóa bỏ các tường xung quanh mô hình. 
2.7.2. Tính toán lực tiếp xúc giữa phần tử và phần tử, giữa phần tử và biên 
Bao gồm 3 bước: 
 Bước 1: tìm kiếm tất cả các phần tử có tiếp xúc với nhau. 
 Bước 2: tính toán lực tiếp xúc lên phần từng phần tử thông qua tiếp xúc. 
Bước 3: tính toán gia tốc, vận tốc, vị trí mới của từng phần tử. 
KiÓm tra tiÕp xóc
(PhÇn tö vµ biªn)
t=
t+
t
TÝnh to¸n lùc tiÕp xóc
TÝnh to¸n vËn tèc, gia tèc
vµ vÞ trÝ cña phÇn tö
Hình 2.23 Trình tự tính toán lực tiếp xúc 
2.7.3. Tích phân phương trình chuyển động 
 Để xác định được lực và chuyển vị của mỗi phần tử theo từng bước thời 
gian cần phải giải phương trình chuyển động [18]. Muốn vậy, phải thiết lập các 
thành phần của véc tơ bao gồm: M , C , F . Quá trình tính toán được chia 
thành nhiều bước thời gian, trong mỗi bước tích toán phải xác định trạng thái 
tiếp xúc của các phần tử trong hệ. Giá trị xuất ra của bước thời gian này lại là 
giá trị đầu vào cho bước thời gian tiếp theo. Do đó, trong mỗi bước tích phân 
các công việc trên đều phải được thực hiện cho đến bước tính toán cuối cùng. 
 76 
Thuật toán phương pháp sai phân trung tâm 
A. Tính toán số liệu ban đầu: 
 1. Nhập số liệu các ma trận:  M và  C 
 2. Xác định các điều kiện ban đầu: 0U và 0U 
 3. Tính gia tốc tại thời điểm ban đầu: 
 
    
 
0 0
0
F C U
U
M

 
 4. Xác định số gia t và tính các hệ số tích phân không đổi: 
 0 2
1
a
t
, 1
1
a
2 t
, 2 0a 2a , 3
2
1
a
a
 5. Tính:    31 0 0 0U U t U a U   
 6. Tính:    0 1M a M a C 
B. Tính toán ở từng bước tích phân: 
 1. Tính toán véc tơ lực tại bước tích phân i: 
          2 2i i i 1i
2 1 1
F(t) F(t) M U M C U
2 tt t 
 2. Giải phương trình với ẩn số là chuyển vị tại bước tích phân i+1: 
  i 1 iM U F(t) 
 3. Xác định vận tốc và gia tốc tại bước tích phân i: 
   1 i 1 i 1iU a U U  
    0 i 1 i i 1iU a U 2 U U  
 4. Xác định vận tốc tại nửa bước tích phân i: 
   1 1 i 1 ii
2
U 2a U U
 và   1 1 i i 1i
2
U 2a U U
 
 77 
2.7.4. Tìm kiếm tiếp xúc 
 Thuật toán tìm kiếm tiếp xúc không nhị phân NBS (No Binary Search) 
được Munjiza và Andrews [16] nghiên cứu để tìm kiếm tiếp xúc giữa các phần 
tử rời rạc có cùng kích thước trong từng bước thời gian. Đối với hệ bao gồm 
nhiều phần tử có kích thước như nhau thì NBS là thuật toán hiệu quả nhất với 
thời gian tính toán ngắn và tiết kiệm được bộ nhớ trong của máy tính. 
Y
X
xmin xmax
ymin
ymax
2
3
4
5
6
7
8
1
Hình 2.24 Mô hình tìm kiếm tiếp xúc 
 Phạm vi nghiên cứu của mô hình nằm trong giới hạn  min maxx x và 
 min maxy y . Thuật toán tìm kiếm tiếp xúc NBS dựa trên việc phân chia miền 
nghiên cứu thành các ô hình vuông giống nhau có kích thước 2R: 
(M,1) (M,2) (M,3) (M,4) (M,5) (M,N) 
(2,1) (2,2) (2,3) (2,4) (2,5)  (2,N) 
(1,1) (1,2) (1,3) (1,4) (1,5)  (1,N) 
Hình 2.25 Phân chia miền tìm kiếm thành các ô hình vuông 
 Số lượng ô vuông là số nguyên được xác định theo 2 phương như sau: 
 max min
x x
ncelx
2R
 và max min
y y
ncely
2R
 (2.133)
 78 
 Tổng số các ô vuông: ncel=ncelx.ncely (2.134)
 Mỗi phần tử được gán với một số nguyên i, i= 1,2,3,..., N (2.135) 
 Số hiệu vị trí các ô vuông (ix,iy): 
 ix=1,2,3,,ncelx và iy=1,2,3,,ncely (2.136) 
 Kiểm tra vị trí của phần tử có tọa độ tâm (x,y) thuộc ô vuông i: 
min
x
x x
i
2R
 và miny
y y
i
2R
 (2.137) 
 Thiết lập địa chỉ cho các phần tử thực hiện theo 2 bước sau đây: 
(Xix-1,Yiy) (Xix,Yiy) 
(Xix-1,Yiy-1) (Xix,Yiy-1) (Xix+1,Yiy-1) 
Hình 2.26 Sơ đồ tìm kiếm tiếp xúc 
 Tìm kiếm tiếp xúc được thực hiện bằng cách kiểm tra địa chỉ của các 
phần tử trong ô vuông đang xét với các phần tử thuộc ô vuông lân cận. Như