Luận án Nghiên cứu chế tạo và tính toán mô phỏng màng mỏng với nền cacbon giống kim cương có tính chống mòn cao và ma sát thấp

 chia khoảng dịch chuyển của mũi kim cương ứng với giai đoạn giảm tải 
cho dịch chuyển lớn nhất. Độ cứng tiếp xúc của màng được xác định bởi đường cong 
gia tải – dịch chuyển trong quá trình đo với lực tăng dần và phân tích bằng phương 
pháp Oliver-Pharr [81] (Hình 4.1, 4.2). 
 Mũi kim cương Berkovich có thể được mô hình hóa bằng một đầu đo hình nón 
với một nửa góc quay Ф=70.30. Chu vi ngoại biên tiếp xúc tạo vết lõm chìm vào vật 
thể được diễn tả bằng mô hình đầu đo cứng có hình học đơn giản [82]. Nếu bỏ qua 
sự tích tụ của vật liệu ở chu vi ngoại biên đối với vật liệu đàn hồi – dẻo, lượng “chìm 
vào” vật liệu màng hs được xác định bằng [81]: 
 푃
 ℎ =∈ (4.1) 
 푠 푆
 43 
trong đó ϵ là hằng số phụ thuộc vào hình dạng của đầu đo (ϵ=0.72 cho hình nón, 0.75 
cho parabol tròn xoay và 1.00 cho dạng phẳng), Pmax tải lớn nhất, S là độ cứng tiếp 
xúc (S=dP/dh). 
 Hình 4.1. Đường cong gia tải – giảm tải [83] 
 Hình 4.2. Dạng hình học của vết lõm khi gia tải và dỡ tải 
 Vì vậy, độ sâu dọc theo lớp tiếp xúc được tạo ra giữa đầu đo và mẫu (hmax-hs) 
là: 
 P
 h = h − ϵ max (4.2) 
 c max S
 Nếu gọi F(d) là “hàm diện tích” miêu tả diện tích được chiếu (hoặc mặt cắt) của 
đầu đo ở khoảng cách d từ bề mặt đến mũi của đo, diện tích tiếp xúc A sẽ là: 
 A = F(hc) (4.3) 
 44 
 Hàm diện tích là hàm hình dạng đầu đo, phải được hiệu chỉnh bằng các phép đo 
độc lập để độ lệch từ hình dạng đầu đo hình nón được đưa vào tính toán thay cho mũi 
Berkovich là nhỏ nhất. Độ lệch này có thể là khá lớn so với mũi kim cương Berkovich 
nếu có sự bo tròn xảy ra trong quá trình mài sắc. Đây là phương pháp cơ bản cho việc 
xác định hàm diện tích tiếp xúc. Pharr và cộng sự [84] đã biểu diễn A là một hàm như 
sau: 
 8 2−n 2 1/2 1/4 1/128
 = ∑n=0 Cn(hC) = C0ℎ + C1h + C2ℎ + C3ℎ + ⋯ + C8ℎ (4.4) 
 Để hoàn thành quá trình hiệu chỉnh hàm diện tích, các hệ số khác trong phương 
trình (4.4) phải được xác định bằng sự điều chỉnh đường cong dữ liệu A. Đối với dữ 
liệu thực nghiệm trong nghiên cứu [85], các hệ số cho hàm diện tích nhận được là 
C0=24.65, C1=202.7, C2=0.03363, C3=0.9318, C4=0.02827, C5= 0.03716, C6=1.763, 
C7=0.04102, và C8= 1.881. 
 Độ cứng của màng được tính bằng công thức: 
 P
 H = max (4.5) 
 A
 Chú ý rằng, việc xác định độ cứng dựa vào diện tích tiếp xúc dưới tác dụng của 
lực, nó có thể lệch so với độ cứng truyền thống được đo từ diện tích của vết lõm nén 
còn dư nếu có sự đàn hồi đáng kể trong suốt quá trình giảm tải. Tuy nhiên, điều này 
chỉ quan trọng với các vật liệu với các giá trị E/H rất nhỏ [86]. 
 Mô đun đàn hồi hiệu dụng (Er) được xác theo mối quan hệ hàm diện tích tiếp 
xúc và độ cứng giảm tải: 
 1 S
 = √ (4.6) 
 훽 2 √ 
 Quan hệ giữa mô đun đàn hồi hiệu dụng, mô đun đàn hồi của đầu đo nano và 
mô đun đàn hồi của mẫu: 
 1 1−v2 1−v2
 = + i (4.7) 
 Er E Ei
trong đó với đầu đâm hình nón thì hệ số hiệu chỉnh β=1000, vật liệu đầu đâm bằng 
kim cương có mô đun đàn hồi E =1040 GPa và hệ số Poisson v = 0,07. Mô đun đàn 
hồi của mẫu (Ei) được xác định khi kết hợp phương trình (4.6) và (4.7). 
 Phương trình (4.7) là công thức rất tổng quát áp dụng cho bất kỳ đầu đo đối 
 45 
xứng trục. Phương trình này không hạn chế đối với dạng hình học đơn giản cụ thể 
nào, kể cả vết lõm tù và có thể áp dụng với tiếp xúc đàn hồi – dẻo [87]. 
4.3. Cơ sở phương pháp phần tử hữu hạn (PTHH) xác định tính chất màng chịu 
tải trọng tĩnh 
4.3.1. Phương trình PTHH 
 Vì có ba ẩn chuyển vị chưa biết u(x,y,z), v(x,y,z) và w(x,y,z) , nên cần 
ba hàm nội suy riêng biệt cho mỗi ẩn. Trên nguyên tắc đó ta có thể sử dụng các hàm 
nội suy khác nhau cho các ẩn khác nhau. Tuy vậy, thực tế người ta sử dụng các hàm 
nội suy giống nhau cho cả ba ẩn chuyển vị như sau [88]: 
 u(x,y,z) = N1u1 + N 2u2 + N 3u3 + 
 v(x,y,z) = N1v1 + N 2v2 + N 3 v3 + 
 w(x,y,z) = N w + N w + N w + 
 1 1 2 2 3 3 (4.8) 
 Ở đây u1, v1, w1, u2,  gọi là các bậc tự do tại nút, Ni (x, y, z) là các hàm nội 
suy tương ứng. Ở dạng ma trận chúng ta có thể viết như sau: 
 u1 
 u N1 0 0 N2 0  v1 
 T
 u(x,y,z)  v  = 0 N1 0 0 N2  . w1  N d
 w 0 0 N 0 0  u 
  1 2 
  
  (4.9) 
 Từ đó có thể tìm được các thành phần biến dạng, ứng suất như sau: 
 u  N N 
 1 0 0 2 0 0 
 x x x 
 v N N 
 0 1 0 0 2 0  
  x  u1 
 y y y 
  v
 y w N1 N2 1 
 0 0 0 0  
 z z z z w1 T
  = = = .  B d
   u v  N N N N    
  xy + 1 1 0 2 2 0  u2 
  y x y x y x v 
 yz 2 
 v w N1 N1 N2 N2 
  zx  + 0 0    
 z y z y z y 
 u w N1 N1 N2 N2
 + 0 0  
 z x z x z x
  (4.10) 
 T
 = C = CB  d  (4.11) 
 46 
trong đó [C] là ma trận các hằng số vật liệu, [B] là ma trận đạo hàm trong quan hệ 
động học giữa chuyển vị - biến dạng, {d} là véc tơ chuyển vị. 
 Năng lượng biến dạng trên phần tử: 
 1 1 T 1 1
 U = T dV = BT d CBT ddV = dT BCBT dV d = dT k dT
 2 2 2 2
 V V V 
 (4.12) 
 Ở đây [k] là ma trận cứng phần tử: 
 [k] = [B][C][B]T dV
 V (4.13) 
 Giả sử rằng các tải trọng tập trung chỉ tác dụng tại các nút, lúc đó công thực 
hiện bởi nó có thể tính như sau: 
 u1 
 v1
 W = F u + F v + F w = F F F F   r T d
 f ( x i i y i i z i i ) x1 y1 z1 x2  w1 f  
 i u 
 2 
   
 (4.14) 
 T
 rf  : chuyển vị của véc tơ tải trọng tác dụng tại nút. 
 T
 rf = Fx1 Fy1 Fz1 Fx2 
 Phiếm hàm thế năng trên phần tử có dạng: 
 1 T T T T
  p (u,v,w) = U − W = d  k  d − rq + rb + rf   d 
 2 (4.15) 
 Từ các điều kiện cần để cực tiểu phiếm hàm thế năng (4.15),   p / di = 0 và 
do các tải trọng tập trung tại nút có thể được lắp ghép một cách trực tiếp vào véc tơ 
tải toàn cục, nên các phương trình cân bằng tĩnh của phần tử như sau: 
 k  d  = rq+ rb
 (4.16) 
 Khi có biến dạng trước  0 do sự thay đổi nhiệt độ hay một lý do nào khác, các 
quan hệ ứng xử (4.11) trở thành: 
  = C  − 
     0 (4.17) 
 47 
 Lưu ý rằng, sự gia tăng nhiệt độ T sẽ dẫn đến sự gia tăng biến dạng đều phụ 
thuộc vào hệ số giãn nở nhiệt của vật liệu, hay nói khác đi sự thay đổi nhiệt độ không 
gây ra biến dạng cắt. Và véc tơ biến dạng trước do sự thay đổi nhiệt độ sẽ là: 
 T
 0 = T T T 0 0 0 (4.18) 
 Kết hợp với các điều kiện ban đầu kèm với ma trận ứng xử [C] là một ma trận 
đối xứng nên biểu thức năng lượng trên phần tử có thể viết lại như sau: 
 1
 U = { }T [C] { } dV − { }T [C] { } dV
 2 0
 V V (4.19) 
 Với lời giải giả định ta có: 
 1 1
U = {d}T [B] [C] [B]T dV{d} − { }T [C] [B]T dV {d} = {d}T [k] {d} − {r }T {d}
 2 0 2 
 V V
 (4.20) 
 T
{r } là chuyển vị của véc tơ tải trọng nút tương đương do có mặt biến dạng ban đầu 
{ 0 } 
 {r } = [B] [C] {0 }dV
 V (4.21) 
4.3.2. Bài toán phẳng của lý thuyết đàn hồi 
 4.3.2.1. Bài toán ứng suất phẳng: 
 Với các vật thể mỏng (giả sử kích thước theo phương z rất nhỏ so với hai phương 
còn lại) chịu tác dụng của các lực trong mặt phẳng xy (Hình 4.3), có thể chấp nhận 
giả thiết [89]. 
 Hình 4.3. Kết cấu tấm cho bài toán ứng suất phẳng 
 48 
 휎 = 휏 = 휏 = 0 (4.22) 
 Lúc này, phương trình quan hệ biến dạng - ứng suất là: 
 1  
 − 0 
  x  E E  x 
  1 
  y  = − 0 .  y 
 E E 
   
 xy  1 xy  (4.23) 
 0 0 
 G 
  + 
  = − x y ;  = 0;  = 0.
 z E yz zx
 Quan hệ ứng suất – biến dạng 
 휎 1 휈 0 휀 
 { 휎 } = [휈 1 0 ] . { 휀 } ⇔ {휎} = [ ]{휀} (4.24) 
 1−휈2 1−휈
 휏 0 0 훾 
 2
 Nếu có sự thay đổi nhiệt độ T , véc tơ biến dạng ban đầu có dạng: 
 T  x +  y
 0 = T T 0 ; z = − + T
 E (4.25) 
   = C  − 0 (4.26) 
 Hình 4.4 thể hiện các ứng suất trong mặt phẳng: 
 Hình 4.4. Nội lực trong một phần tử tấm ứng suất phẳng 
4.3.2.2. Bài toán biến dạng phẳng: 
 Với các vật thể dài (giả sử kích thước theo phương z rất lớn so với hai phương 
 49 
còn lại) chịu tác dụng của các lực trong mặt phẳng xy, người ta có thể chấp nhận giả 
thiết 
  =  =  = 0
 z yz zx (4.27) 
 Từ đó, ta có quan hệ ứng suất – biến dạng 
 1 − 휈 휈 0
 휎 휀 
 휈 1 − 휈 0
 { 휎 } = [ ] . { 휀 } 
 (1 + 휈)(1 − 2휈) 1 − 2휈
 휏 0 0 훾 
 2
 휀 +휀
 휎 = 휈 = 0; 휏 = 0; 휏 = 0 (4.28) 
 (1+휈)(1−2휈) 
 4.3.2.3. Các phương trình cân bằng: 
 Đối với bài toán phẳng, việc giả sử chiều dày h theo phương z là hằng số dẫn 
đến các tích phân thể tích sẽ trở thành các tích phân trên diện tích phần tử, các tích 
phân mặt sẽ trở thành các tích phân trên biên phần tử. Tải trọng tác dụng chỉ có các 
thành phần trong mặt phẳng xy. Do đó các phương trình phần tử cho bài toán phẳng 
có thể thiết lập như sau: 
 Quan hệ biến dạng – chuyển vị 
 u  N N 
 1 0 2 0 
 u1 
  x  x x x 
 v N N v 
  =  = = 0 1 0 2  . 1  [B]T d
  y    
 y y y u2 
  xy  u v N N N N 
 1 1 2 2  
 +  
 y x  y x y x (4.29) 
 Quan hệ ứng suất – biến dạng 
  = C  − 
     0 (4.30) 
 Cho bài toán ứng suất phẳng: 
   
  x 1  0  x T
 E 
  =  1 0 .  − T
 y  2 y  
 1 − 1 − (4.31) 
  xy 0 0  xy 0 
  2  
 Cho bài toán biến dạng phẳng: 
 50 
   1 −  0   (1 + ) T
 x x 
 E 
  y  =  1 − 0 .  y  − (1 + ) T (4.34) 
 (1 + )(1 − 2 ) 1 − 2 
  z  0 0  xy 0 
 2  
 Phương trình cân bằng của phần tử 
 k  d = r + r + r 
 q  b  (4.35) 
trong đó, ma trận cứng phần tử được xác định 
 k  = h [B][C][B]T dA
 A (4.36) 
 Véc tơ tải trọng nút do tải trọng phân bố 
 푞 
 { 푞} = ℎ ∫ [ ] { } (4.37) 
 푞 
 [N]: là ma trận hàm dạng. 
 Véc tơ tải trọng do tải thể tích 
 { } = ℎ ∬ [ ] { } (4.38) 
 Trong luận án, tác giả sử dụng phần tử tứ giác 4 điểm nút để mô hình hóa kết 
cấu màng mỏng. Phương trình hàm dạng [N] trong hệ tọa độ tự nhiên của phần tử tứ 
giác 4 điểm nút [90]: 
 1
 N = (1−  )(1−)
 1 4
 1
 N = (1+  )(1−)
 2 4
 1
 N = (1+  )(1+)
 3 4
 1
 N 4 = (1−  )(1+)
 4 (4.39) 
 Mối liên hệ hình học và nội suy chuyển vị được biểu diễn bởi: 
 1 1 1 1 1
 푒(휉, 휂)
 1 
 1 2 3 4 푒(휉, 휂)
 2 
 = 1 2 3 4 푒 (4.40) 
 3 (휉, 휂) 
 1 2 3 4 
 [ 푒(휉, 휂)]
 [ ] [ 1 2 3 4] 4
 51 
4.3.3. Tích phân Gauss 
 Tích phân Gauss được viết dạng tổng quát: 
 푒 +1 +1[ 푒] 푒 | |
 = ∫−1 ∫−1 휌 푒푡 퐽 휉 휂 (4.45) 
 Trong đó ma trận hàm dạng phần tử Ne có dạng: 
 푒 푒( ) 푒( ) 푒 ( )
 (휉, 휂) = [ 1 휉, 휂 2 휉, 휂  푛푛 휉, 휂 ] (4.46) 
 Công thức Gauss cho trường hợp hai chiều, với điều kiện miền lấy tích phân là 
một miền vuông có cạnh 2 đơn vị, gốc tọa độ đặt tại tâm miền như Hình 4.5. Lúc đó 
chúng ta khảo sát trên tích phân sau: 
 1 1
 I = f (x,y) dx dy
 −1 −1 (4.47) 
 Hình 4.5. Miền lấy tích phân cho công thức Gauss hai chiều 
 Xét dải tô đen theo phương đứng như Hình 4.5, với mỗi x bất kỳ cho trước, tích 
phân theo phương y có thể tính được bằng cách sử dụng công thức Gauss một chiều 
n điểm. Và tương tự sử dụng công thức Gauss với m điểm theo phương x. Cuối cùng 
tích phân trên có thể viết lại như sau: 
 1 n m n
I w j f (x,y j ) dx wi w j f (xi ,y j )  wi w j f (xi ,y j )
 −1  j=1 i=1  j=1 
 i, j
 (4.48) 
 Việc tìm tọa độ và trọng số tương ứng của các điểm Gauss cũng tương tự như 
đối với trường hợp một chiều, Bảng 4.1 cho ta giá trị của một vài trường hợp cụ thể. 
 52 
 Bảng 4.1. Trọng số tương ứng của các điểm Gauss [88] 
 x j y j wi w j
1x1 1 0 0 4.0000000000000 
 1 - 0.57735026918963 - 0.57735026918963 1.000000000000 
 2 - 0.57735026918963 0.57735026918963 1.000000000000 
2x2 
 3 0.57735026918963 - 0.57735026918963 1.000000000000 
 4 0.57735026918963 0.57735026918963 1.000000000000 
 1 - 0.77459666924148 - 0.77459666924148 0.30864197530864 
 2 - 0.77459666924148 0 0.49382716049383 
 3 - 0.77459666924148 0.77459666924148 0.30864197530864 
 4 0 - 0.77459666924148 0.49382716049383 
3x3 5 0 0 0.7901234567901 
 6 0 0.77459666924148 0.49382716049383 
 7 0.77459666924148 - 0.77459666924148 0.30864197530864 
 8 0.77459666924148 0 0.49382716049383 
 9 0.77459666924148 0.77459666924148 0.30864197530864 
 1 - 0.86113631159405 - 0.86113631159405 0.12100299328560 
 2 - 0.86113631159405 - 0.33998104358486 0.22685185185185 
 3 - 0.86113631159405 0.33998104358486 0.22685185185185 
 4 - 0.86113631159405 0.86113631159405 0.12100299328560 
4x4 5 - 0.33998104358486 - 0.86113631159405 0.22685185185185 
 6 - 0.33998104358486 - 0.33998104358486 0.42529330301069 
 7 - 0.33998104358486 0.33998104358486 0.42529330301069 
 8 - 0.33998104358486 0.86113631159405 0.22685185185185 
 9 0.33998104358486 - 0.86113631159405 0.22685185185185 
 53 
 w w
 10 0.33998104358486x j - 0.33998104358486y j 0.42529330301069i j 
 11 0.33998104358486 0.33998104358486 0.42529330301069 
 12 0.33998104358486 0.86113631159405 0.22685185185185 
 13 0.86113631159405 - 0.86113631159405 0.12100299328560 
 14 0.86113631159405 - 0.33998104358486 0.22685185185185 
 15 0.86113631159405 0.33998104358486 0.22685185185185 
 16 0.86113631159405 0.86113631159405 0.12100299328560 
4.4. Phương trình ứng xử của màng mỏng chịu tải bởi đầu đâm nano và phương 
pháp giải 
 Như trình bày ở mục 4.1, vết lõm là một phương pháp đơn giản và thường được 
áp dụng để kiểm tra tính chất cơ của các vật liệu. Trong kiểm tra vết lõm, vật liệu 
quan tâm chưa được biết các tính chất cơ học (như mô đun đàn hồi và độ cứng) được 
tiếp xúc với vật liệu đã biết được các tính chất cơ học. Kỹ thuật bắt nguồn từ thang 
độ cứng Mohs được giới thiệu năm 1812. Đối với việc kiểm tra vết lõm truyền thống, 
chiều sâu thâm nhập được đo ở mức độ μm hay mm. Gần đây hơn, vết lõm nano (vết 
lõm mà có chiều sâu thâm nhập được đo ở mức độ nano) đã được thiết lập như một 
công cụ cơ bản cho việc nghiên cứu, kiểm tra tính chất cơ học của vật liệu thể tích 
nhỏ [89 - 92]. Kỹ thuật vết lõm nano có thể được dùng không chỉ để tính toán độ 
cứng, mô đun đàn hồi mà còn tính các chỉ số biến dạng – sự hóa cứng, độ bền chống 
nứt (cho vật liệu giòn) và các tính chất đàn hồi. Trong nghiên cứu này, trường chuyển 
vị, biến dạng, ứng suất phân bố trong màng sẽ được tính toán và phân tích. 
 Việc xây dựng mô hình toán màng mỏng chịu tải bởi đầu đâm nano (Hình 4.6) 
với giả thiết màng mỏng nanocomposite với các hạt tinh thể TiC phân tán trong nền 
vô định hình DLC phủ lên đế Si [100] có đường kính 100 mm, dày 450 μm. Màng 
được phủ lên đế với độ dày 1,3 μm. Giả thiết liên kết giữa màng và đế là lý tưởng, 
quá trình chịu tải của màng không bị ảnh hưởng bởi ma sát, nhiệt độ, các hạt tinh thể 
TiC phân tán đồng đều trong nền DLC. Đầu đâm nano là mũi kim cương Berkovich 
 54 
có mô đun đàn hồi 1040 GPa, hệ số poisson 0,07, đâm xuống màng theo phương 
thẳng đứng và tăng dần sau mỗi bước. Tuy nhiên, để thuận lợi trong tính toán, mũi 
Berkovich được qui đổi về dạng đối xứng trục hình nón tương đương có nửa góc 
70.30 [93]. 
 Do quan hệ đường cong gia tải – giảm tải như hình hình 4.1, nên mô đun đàn 
hồi vật liệu màng thay đổi theo tải trọng tác dụng, do đó ma trận [k] trong phương 
trình (4.35) phụ thuộc chuyển vị - phương trình phi tuyến vật lý. Điều này dẫn đến 
việc giải phương trình này được thực hiện trên từng bước gia tải, tương ứng với các 
số gia chuyển vị. Lúc này phương trình (4.35) được viết cho số gia như sau: 
 [퐾]{∆ }𝑖 = {∆ }𝑖 (4.49) 
 [K] - ma trận độ cứng tại bước gia tải { f}i, 
 { u}i - số gia véc tơ chuyển vị của nút tại bước i, 
 { f}i - số gia véc tơ lực 
 i chỉ số thứ tự bước tính. 
 Hình 4.6. Mô hình màng chịu tải 
 Để có được đáp án cho giới hạn biên, các thay đổi của giới hạn biên được 
liên tục cập nhật theo các bước tiến tại công thức (4.49). Kết quả cuối cùng là tổ 
hợp của tất cả kết quả của các bước trước đó. Tính phi tuyến của vật liệu thể hiện 
qua ma trận độ cứng chung, thông số này phụ thuộc vào ứng suất hiện tại và độ 
 55 
biến dạng và thay đổi theo mỗi bước tiến, do đó ma trận độ cứng của màng được 
cập nhật sau mỗi bước tính. Sơ đồ của phương pháp được thể hiện như hình 4.8. 
 (a) độ cứng không đổi (b) độ cứng tiếp tuyến 
 Hình 4.7. Các phương pháp phi tuyến tính [94] 
 Từ (4.49), ta có: 
 {∆ }𝑖 = [퐾]−1{∆ }𝑖 (4.50) 
 {∆휀}𝑖 = [B] {∆ }𝑖 (4.51) 
trong đó [B] là ma trận đạo hàm trong quan hệ chuyển vị - biến dạng. 
 Nếu vật liệu biến dạng dẻo, biến dạng bao gồm đàn hồi và dẻo: 
 {∆휀}𝑖 = {∆휀푒}𝑖 + {∆휀 }𝑖 (4.52) 
 Lúc này ứng suất trong màng được xác định: 
 {∆휎}𝑖 = [ ]{∆휀푒}𝑖 (4.53) 
 Và do đó: 
 {휎}𝑖 = {휎}𝑖−1 + {∆휎}𝑖 (4.54) 
 Phương trình giới hạn chảy được xác định như sau [95]: 
 휎 휎 휎 휎
 퐹 = ( 1+ 3) 푠𝑖푛φ − ( 1− 3) − . 표푠φ = 0 (4.55) 
 2 2
trong đó c là độ bền trượt, φ là góc ma sát của sai số Mohr Coulomb 
 Một trong những đặc điểm của bài toán luận án giải quyết là nếu ứng suất vẫn 
dưới trạng thái giới hạn này (biến dạng đàn hồi), do quan hệ tải trọng - chuyển vị là 
đường cong, nên việc giải lặp theo (4.49) vẫn phải được thực hiện. 
 Trong trường hợp ứng suất vượt quá giới hạn (phương trình (4.55), F>0), ứng 
 56 
suất dư cần được cân bằng lại để đạt tới giá trị gần bằng với giới hạn (F=0). Khi đó 
vật liệu bị chảy dẻo và giới hạn chảy được tính theo Vermeer [96]: 
 휕 𝑖−1 휕 𝑖
 {∆휀 } = ∆휆 [(1 − 휔) ( ) + 휔 ( ) ] (4.56) 
 휕휎 휕휎
trong đó, ∆휆 là số gia biến dạng dẻo và là tham số phụ thuộc vào thời gian sử dụng. 
 Với 휔 =0 thì phương trình được gọi là hàm tường minh và 휔=1 thì gọi là hàm 
ẩn. Sử dụng hàm tường minh với ứng suất lấy từ bước trước 휔 =0, do đó công thức 
(4.56) được viết thành: 
 휕 𝑖−1
 {∆휀 } = ∆휆 ( ) (4.57) 
 휕휎
 Biến dạng dẻo được tính dựa trên ứng suất “dư” tại các nút và ứng suất này được 
cộng vào bước tiếp theo theo luật lũy kế. Lực tạo ra được tính theo phương pháp biến 
dạng dẻo Visco [97]. 
 Trong quá trình tính sử dụng các phương trình từ (4.49) tới (4.57) được lặp đi 
lặp lại cho tới khi ứng suất gần bằng với giới hạn chảy bề mặt trong khoảng sai số 
cho phép nhỏ hơn 0,01 thì bài toán được xem như hội tụ. Các điều kiện hội tụ gồm: 
 + Bước chuyển vị từ phương trình (4.50) gần bằng với bước liền kề. 
 + Lượng biến dạng quá nhỏ. 
 + Ứng suất tăng từ phương trình (4.54) gần bằng với bước liền kề. 
 Nếu sử dụng điều kiện chuyển vị để tính toán sai số thì giá trị chuyển vị của 
bước tiến thứ (i) và bước tiến thứ (i + 1) phải nhỏ hơn sai số cho phép. Với bài toán 
luận án tác giả sử dụng điều kiện hội tụ theo ứng suất. 
 Lưu đồ thuật giải bài toán được thể hiện như trên Hình 4.8. Việc cụ thể hóa thuật 
giải được tác giả lập trình trong môi trường Matlab. Chương trình tính chạy trên máy 
tính với hệ điều hành Windows, có khả năng phân tích tĩnh học phi tuyến của màng 
mỏng. Cấu trúc của chương trình gồm các mô đun chính: 
 Mô đun nhập số liệu: data.m; 
 Mô đun xây dựng mô hình hình học: model.m; 
 Mô đun phân tích tĩnh: static_analysis.m; 
 Mô đun xuất số liệu: results.m. 
 57 
 Bắt đầu 
 Nhập số liệu 
 Thiết lập thông số ban đầu và điều kiện biên 
 Thiết lập [C], [K], {σ0}, {ε0} 
 Bắt đầu bước tính j 
 Đọc bước tăng tải {Δf}i Tái lập bước tiến chuyển vị 
 i 
 Tái lập {Δr}
 Bắt đầu bước 
 Cập nhật tải {Δf}i 
 i Cập nhật diện tích mặt biến dạng 
 Bước tiến chuyển vị {Δu} Cập nhật lực nén bề mặt P 
 Cập nhật bước tiến chuyển vị 
 dui 
 Tổng biến dạng {Δε}i Cập nhật ứng suất {σ}j 
 Tổng ứng suất {Δσ}i Tính phương trình biến dạng dẻo F 
 Tới điểm Sai 
 F > 0 
Gauss tiếp Đúng 
 theo 
 Bước tăng biến dạng dẻo {δεp}i 
 Cập nhật biến dạng dẻo {Δεp}i 
 Tính toán tải tại điểm {Δr}i 
 Kết thúc 
 Thiết lập tất cả các điểm {Δr}i 
điểm Gauss 
 Tính toán sai số e 
 Sai 
 e< giớ i hạn 
 Đúng 
 j j
 Cập nhật chuyển vị {u} , Cập nhật ứng suất {σ} , 
 j 
 Cập nhật biến dạng {ε} , Cập nhật độ sâu đâm h 
 Sai 
 h > chi ều sâu 
 Đúng 
 Hồi phục 
 Thiết lập chuẩn năng lượng biến dạng U, sao lưu dữ liệu, xuất biểu đồ 
 Kết thúc 
 Hình 4.8. Lưu đồ thuật giải 
 58 
4.5. Mô phỏng ứng xử của màng mỏng chịu tải bởi đầu đâm nano 
 Bảng 4.2. Thông của màng và đầu đâm nano, đế Si[100] 
 Mô đun đàn hồi ban 
 Kích thước Poisson 
 đầu (GPa) 
 Đầu đâm nano Hình nón 1040 0,07 
 Màng 
 Bề dày 1,3 μm 241 0,23 
 nanocomposite 
 D=100 mm, 
 Đế 185 0,28