Khảo sát ảnh hưởng của sự phân bố lỗ rỗng tới sự khởi tạo và phát triển vết nứt bằng phương pháp trường pha

100 TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
KHẢO SÁT ẢNH HƯỞNG CỦA SỰ PHÂN BỐ LỖ RỖNG 
TỚI SỰ KHỞI TẠO VÀ PHÁT TRIỂN VẾT NỨT 
BẰNG PHƯƠNG PHÁP TRƯỜNG PHA
Nguyễn Thị Hải Như1*, Trần Anh Bình2 
Tóm tắt: Các loại vật liệu tựa dòn thường được xem là đồng nhất ở tỷ lệ lớn, nhưng gồm nhiều pha ở tỉ lệ 
nhỏ trong đó có các lỗ rỗng và các vết nứt siêu nhỏ. Các vết nứt và lỗ rỗng này có vai trò quan trọng trong 
việc hình thành vết nứt thực ở tỷ lệ lớn hơn. Sự phân bố của các vết nứt và các pha rỗng ở tỷ lệ nhỏ có ảnh 
hưởng khác nhau đến sự khởi tạo và lan truyền vết nứt trong vật liệu. Nhiều mô hình đã được xây dựng để 
nghiên cứu sự lan truyền của vết nứt. Gần đây, các mô hình trường pha nhận được rất nhiều sự quan tâm 
nghiên cứu. Bài báo này xây dựng công thức phần tử hữu hạn cho một mô hình trường pha, sử dụng mô 
hình để khảo sát ảnh hưởng của phân bố lỗ rỗng tới sự phá hoại của một mẫu vật liệu chịu kéo.
Từ khóa: Khởi tạo vết nứt; phát triển vết nứt; phương pháp trường pha.
Study the impact of voids distribution on the initiation and propagation of cracks by a phase field method
Abstract: Quasi brittle materials are considered to be homogeneous at macroscale but composed of phases 
at microscale including void inclusions and microcracks. These defects play an important role in the forming 
a real crack at macroscale. The distribution of microcracks and voids affects differently on the initiation and 
propagation of cracks. A lot of works have been proposed to study the initiation and propagation of the crack 
in literature. Recently phase field models are of interest by many researchers who work on the damage of 
materials because of the ability to track the crack naturally by solving the minimization problem of energy. 
This work will derive discretization formulas for a phase field model and ultilize this model to study the impact 
of the void distribution on the fracture process of a specimen in a traction test.
Keywords: Crack initiation; crack propagation; microvoid; phase field method.
Nhận ngày 22/3/2017; sửa xong 20/6/2017; chấp nhận đăng 26/9/2017 
Received: March 22th, 2017; revised: June 20th, 2017; accepted: September 26th, 2017
1. Giới thiệu
Các loại vật liệu xây dựng như bê tông, 
đá, ceramic, vữa, thạch cao thường tồn tại các 
pha rỗng và các vết nứt siêu nhỏ không thể thấy 
được bằng mắt thường. Các hư hại này được 
hình thành một cách tự nhiên trong quá trình hình 
thành nên vật liệu (ví dụ: Hình 1 thể hiện các vết 
nứt ở tỷ lệ micro của vật liệu bê tông nhẹ, hình 
thành do quá trình thủy hóa, đông cứng của bê 
tông). Dưới tác dụng của tải trọng, nhiệt độ, các 
hư hại ở các vị trí khác nhau lan truyền, rẽ nhánh 
và kết hợp với nhau tạo thành các vết nứt thực sự 
có thể nhìn thấy bằng mắt thường và dẫn đến sự 
phá hoại của toàn kết cấu.
Trong nhiều thập niên vừa qua, nghiên cứu về sự hư hại của vật liệu và kết cấu nhận được sự quan 
tâm của nhiều tác giả. Khi nghiên cứu vết nứt bằng phương pháp phần tử hữu hạn truyền thống (FEM), tại 
1 ThS, Khoa Công nghệ thông tin, Trường Đại học Xây dựng.
2 TS, Khoa Công nghệ thông tin, Trường Đại học Xây dựng.
* Tác giả chính. E-mail: nhunth@nuce.edu.vn.
Hình 1. a) Vi cấu trúc trong bê tông nhẹ 
(Amine Bouterf-phD thesis 2014) và b) các vết nứt ở 
tỉ lệ micro (Passmore et al.1965)
101TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
vị trí bị phá hoại, các nút phải được tách rời dẫn đến việc chia lưới phải được cập nhật và đường đi của vết 
nứt được quy định bởi cạnh phần tử [1]. Việc này có thể giải quyết bằng các thuật toán chia lại lưới trong 
[2,3] hoặc sử dụng phương pháp XFEM [4-7] trong đó sử dụng các hàm bình đồ (hàm level-set) để biểu diễn 
sự không liên tục trong phần tử. Các phương pháp này đều có những thành công nhất định trong mô hình 
vết nứt 2 chiều nhưng gặp khó khăn khi vết nứt rẽ nhánh trong mô hình 3 chiều. Để giải quyết khó khăn này, 
một nhóm các phương pháp đã được đề xuất trong đó vết nứt được mô tả là một đại lượng có chiều dày 
hữu hạn. Trong số các phương pháp này trường pha nhận được nhiều sự quan tâm của các nhà nghiên cứu 
hiện nay. Trong các mô hình này, vết nứt được biểu diễn thông qua một trường vô hướng có giá trị từ 0 đến 
1. Cách biểu diễn này tương đồng với cách biểu diễn biến hư hại trong nghiên cứu về hư hại trong cơ học 
môi trường liên tục về hư hại (Continuum Damage Mechanics). Bài toán về phá hoại được giải quyết dựa 
trên nguyên lý biến phân không cần đến các tiêu chuẩn phá hoại cũng như các thuật toán chia lại lưới và 
làm cho việc mô hình vết nứt trong không gian 3 chiều cũng trở nên dễ dàng hơn. Với lợi thế này, rất nhiều 
mô hình trường pha được phát triển để nghiên cứu các bài toán hư hại khác nhau: phá hoại kiểu I [8], kiểu 
III [9], hay kiểu I và II [10]. Trong lĩnh vực cơ học, các phương pháp trường pha xây dựng trên nguyên lý biến 
phân cho vật liệu dòn trong đó sử dụng một hàm trơn, khả vi để mô tả bề mặt nứt đã được đề xuất trong [11-
13]. Loại mô hình này được mở rộng áp dụng cho vết nứt cố kết [17,18] và các bài toán động học [14-16]. 
Bài báo này xây dựng các công thức phần tử hữu hạn để tìm trường hư hại sử dụng mô hình trường 
pha (phase field) [19,20] theo nguyên lý cực tiểu hóa năng lượng ở mục 2. Trong mục 3 trình bày kết quả 
tính toán theo mô hình trong bài báo và so sánh với một kết quả đã được công bố. Trong mục 4, nhóm tác 
giả tiến hành khảo sát ảnh hưởng của sự phân bố các pha rỗng tròn xung quanh vết nứt của một mẫu vật 
liệu chịu kéo. Các kết luận từ việc tính toán khảo sát ở trên được đưa ra trong mục 5.
2. Phương pháp trường pha 
2.1 Công thức năng lượng của mô hình
Xét một miền mở ký hiệu là Ω và biên của nó là dΩ trong không gian D chiều với n là vec tơ pháp 
tuyến hướng ra ngoài như mô tả trong Hình 2. Miền này chứa một vùng nứt Γ trong không gian D - 1 chiều. 
Hình 2. Miền Ω và chứa 
miền nứt Γ
Hình 3. Mô tả giá trị của d trong mô hình vết nứt thực (sharp crack) 
và vết nứt mờ (smear crack) trong bài toán 1 chiều
Francfort và Marigo [11] đề xuất công thức năng lượng cho bài toán mô hình vết nứt như sau:
 (1)
Trong mô hình chuẩn hóa của lý thuyết của trường pha, hàm năng lượng E(u,Γ) được biểu diễn dưới 
dạng sau:
 (2)
Trong (1) và (2): Wu là hàm mật độ năng lượng biến dạng đàn hồi; u là trường chuyển vị; 
là ten-sơ biến dạng, gc là tỷ lệ giải phóng năng lượng tới hạn; Eu đại diện cho năng lượng 
biến dạng đàn hồi, Es là năng lượng bề mặt hay năng lượng phá hủy. Lưu ý trong công thức (1), Es là một 
đại lượng tích phân trên bề mặt nứt, còn trong công thức (2), Es là một đại lượng được tính tích phân trên 
toàn miền của hàm mật độ trường nứt γ(d,∇d). Hàm mật độ được chọn theo [19] như sau:
 (3)
102 TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
Trong mô hình đang xét, d = 1 biểu thị vật liệu ở trạng thái hoàn toàn bị phá hủy tại vị trí vết nứt; 
d = 0 biểu thị trạng thái nguyên vẹn của vật liệu. Trên Hình 3, với bài toán 1 chiều, giá trị của d là một miền 
liên tục d, trong đó d = 1 tại vị trí vết nứt và liên tục giảm dần khi ra ngoài miền hư hại và tiến đến 0. Độ lớn 
của miền phụ thuộc vào biến l trong công thức (3). Tham số l ban đầu được giới thiệu là một tham số của 
mô hình và nên được chọn càng bé càng tốt tuy nhiên thông số này có mối quan hệ với các tính chất khác 
của vật liệu. Do đó, l nên được xem là một tính chất của vật liệu và được xác định dựa theo thí nghiệm.
Công thức (2) có thể được viết lại dưới dạng: trong đó W = W(ε(u),d) + gcγ(d) gọi là 
năng lượng tự do.
Để mô tả sự hư hại của vật liệu, hàm năng lượng biến dạng trong các mô hình trường pha thường 
được viết kèm theo một hàm suy thoái g(d). Dựa trên giả thiết rằng chỉ có biến dạng kéo gây ra vết nứt, hàm 
suy thoái được đưa vào công thức gắn với phần dương của năng lượng biến dạng đàn hồi. Theo đó, công 
thức năng lượng cho mô hình có dạng:
 (4)
Trong công thức (4): k là một hệ số rất nhỏ đế tránh biến dạng kì dị khi phần tử bị phá hoại hoàn toàn 
(d = 1 trên toàn miền); hàm g(d) được chọn đơn giản nhất có dạng g(d) = (1 – d)2 thỏa mãn các điều kiện: 
g(d=0) = 1; và g(d=1) = 0 tương ứng giới hạn 2 trạng thái ban đầu và trạng đã bị hư hại, g'(d=1) = 0 nhằm đảm bảo 
rằng khi hư hại hoàn toàn xảy ra, năng lượng của mô hình hội tụ đến một giá trị hữu hạn và trường hư hại 
không thể tiếp tục tiến triển.
Theo [19,20], năng lượng dương và năng lượng âm được tính dựa trên phần kéo và phần nén của 
ten xơ biến dạng. Các đại lượng trong (4) được tính như sau: 
 (5)
Với
 (6)
Trong các công thức trên, εi và ni là các trị riêng và vec tơ riêng của ten sơ biến dạng ε; 
và .
2.2 Phương trình biến phân
Để giải bài toán cực tiểu năng lượng theo công thức (2), hai giải thuật chính được đề nghị trong các 
nghiên cứu về mô hình bao gồm giải đồng thời tìm giá trị của trường u và d (monothilic scheme) và giải lần 
lượt d và u (staggered scheme) [19,20]. Trong bài báo này, giải thuật thứ 2 được chọn do tính đơn giản và 
cho kết quả hợp lý. 
Công thức dạng yếu của bài toán tìm trường d (phase field) có thể được xây dựng dựa trên sự nhất 
quán với các nguyên lý nhiệt động học (xem [20]). Công thức tích phân để tìm d cũng có thể xây dựng dựa 
theo nguyên lý biến phân trong đó d là nghiệm của bài toán cực tiểu hóa năng lượng theo biến d hay đạo 
hàm theo hướng δd của E bằng 0. Ta có thể viết:
 (7)
Thay (3) vào (7), bài toán tìm trường d được viết lại:
hay (8)
Khai triển (8), ta có:
103TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
 (9)
hay
 (10)
Với giả thiết vết nứt không thể liền lại, Miehe [19] giới thiệu hàm H thay thế cho hàm Ψ+ trong công 
thức (10).
 (11)
Thay (11) vào (10) ta được công thức dạng yếu tìm trường d:
 (12)
Công thức dạng tích phân của bài toán tìm chuyển vị được viết tương tự như trong phương pháp 
phần tử hữu hạn truyền thống:
 (13)
trong đó: Su = {u|u(x) = u trên ∂Ω} là không gian các hàm u thỏa mãn điều kiện biên; Fext = ∫ΩfudΩ + 
∫∂ΩFF̅udΩ), f là lực thể tích; F̅ là lực tác dụng trên biên ∂ΩF và u là chuyển vị cưỡng bức đặt trên biên ∂Ωu với 
điều kiện ∂Ω = ∂Ωu U ∂ΩF. 
Khi đó, phương trình biến phân của bài toán đàn hồi có vết nứt theo phương pháp trường pha được 
viết dưới dạng:
 (14)
trong đó: 
 (15)
Trường biến dạng được tách thành 2 thành phần theo công thức xấp xỉ như sau:
 (16)
trong đó:
 (17)
2.3 Rời rạc hóa bằng phương pháp phần tử hữu hạn
2.3.1 Công thức ma trận giải bài toán tìm trường d
Áp dụng các công thức gần đúng của phương pháp phần tử hữu hạn cho biến d và δd:
 (18)
trong đó: Nd(x) và Bd(x) lần lượt là ma trận hàm dạng và ma trận đạo hàm của hàm dạng. Sau khi biến đổi 
(12) ta thu được hệ phương trình tuyến tính:
Kdd = Fd (19)
trong đó, ma trận độ cứng Kd và vec tơ Fd lực được xác định bởi:
 (20)
104 TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
2.2.2 Công thức ma trận giải bài toán tìm trường chuyển vị
Áp dụng các công thức gần đúng của phương pháp phần tử hữu hạn, ta nhận được:
 (21)
trong đó:
 (22)
Trong công thức (20): λ và μ là các hằng số l’amé của vật liệu, [1] = [1 1 0]. N và B lần lượt là ma trận 
hàm dạng và ma trận đạo hàm của hàm dạng của phép xấp xỉ chuyển vị u. Pn+ và Pn- được giới thiệu như 
trong (17); R+ và R- là các toán tử xấp xỉ thỏa mãn điều kiện sau:
 (23)
hoặc có thể viết: 
Tại mỗi bước chuyển vị, trường d được tìm theo công thức (19) và được sử dụng để giải tìm trường 
chuyển vị u theo công thức (20).
3. Ví dụ tính toán: tấm chịu kéo
Ví dụ này nhằm mục đích so sánh kết quả chương trình phần tử hữu hạn đã lập theo mô hình trình 
bày trong mục 2 với kết quả đã được công bố trong [19].
Xét một mẫu biến dạng phẳng có kích thước L × H = 1mm × 1mm chứa vết nứt ban đầu có chiều dài 
0.5mm tại vị trí Y = 0.5mm. Các hằng số vật liệu λ = 121.15KN/mm2, μ = 80.77 KN/ mm2, gc = 0.0027 KN/
mm. Tấm được ngàm tại đáy như Hình 4.a. Tiến hành áp chuyển vị tăng dần tại vị trí Y = 1mm: ∆u = 10-5 mm 
cho 500 bước chuyển vị đầu tiên và ∆u = 10-6 mm cho các chuyển vị tiếp theo. Hệ số l = 0.015mm, hệ số 
k = 10-6. Lưới chia phần tử hữu hạn được thể hiện trong Hình 4.b trong đó sử dụng 15292 phần tử tam giác, 
kích thước phần tử tại vùng dự đoán hư hại xấp xỉ 0.001mm. 
Hình 4.c thể hiện quan hệ lực chuyển vị được tính toán theo phương pháp trường pha, so sánh với 
kết quả trong [19]. Đồ thị cho thấy kết quả của nhóm tác giả là đáng tin cậy.
Hình 4. a) Sơ đồ hình học; b) Lưới chia phần tử hữu hạn; c) Kết quả tính sử dụng mô hình trường pha trong 
mục 2 (đường liền) được so sánh với kết quả trong [19] (đường nét đứt)
4. Ảnh hưởng của sự phân bố của lỗ rỗng đến sự hình thành và phát triển của vết nứt
Trong phần này, ảnh hưởng của sự phân bố lỗ rỗng đến sự hình thành và phát triển của vết nứt 
được khảo sát bằng phương pháp trường pha qua 4 mẫu biến dạng phẳng có kích thước L × H = 1mm × 
1mm. Các mẫu đều có 11 lỗ tròn nhưng phân bố khác nhau và một vết nứt ban đầu song song với trục X. 
Bán kính các của lỗ r = 0.02mm phân bố tạo thành dạng tam giác đều như trên Hình 5.a, vết nứt ban đầu 
có chiều dài 0.5 mm tại vị trí Y = 0.5mm. Các hằng số vật liệu được lấy như ví dụ trong mục 3, tham số 
l = 0.01mm. Ứng xử của các mẫu được khảo sát dựa trên quan hệ lực - chuyển vị và hình dạng vết nứt 
105TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
(phân bố của trường d) của các thí nghiệm kéo bằng cách áp chuyển vị tăng dần theo phương Y tại vị trí 
Y = 1mm (chuyển vị theo phương X, Ux = 0). Bước nhảy chuyển vị ∆u = 10-4 mm và 120 bước chuyển vị 
được thực hiện trong mỗi mẫu. 
Hình 5. Cột a: Các mẫu thử; Cột b: Phân bố trường hư hại d (vị trí các vết nứt) tại điểm cuối của 
quá trình gia tải; Cột c: Mối quan hệ lực - chuyển vị cưỡng bức
Xấp xỉ 80,000 phần tử tam giác tuyến tính được sử dụng để rời rạc hóa phần tử. Việc chia lưới được 
thực hiện bằng phần mềm ABAQUS 6.14. Việc tính toán được thực hiện bằng phần mềm MATLAB 2009. 
Trên Hình 5.b thể hiện trường phá hoại d và Hình 5.c thể hiện quan hệ lực chuyển vị F̅ – u̅ tương ứng với 
các mẫu.
Quan sát hai mẫu m1 và m2 (Hình 5): Các lỗ được bố trí phía trước vết nứt ban đầu. Tọa độ theo 
phương X ở hai mẫu hoàn toàn giống nhau nhưng khác nhau theo phương Y. Khoảng cách giữa tâm các 
lỗ tương ứng là 0.08mm và d = 0.1386 mm. Ở mẫu m1 các lỗ rỗng được phân bố gần nhau, trường phá 
hoại xuất hiện tại tất cả các mép lỗ và vết nứt tìm đi theo hình dích dắc nối các lỗ với nhau. Ở mẫu m2 các 
106 TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
lỗ rỗng được phân bố xa nhau với khoảng cách 
giữa tâm các lỗ, vết nứt không đi qua tất cả các lỗ 
mà đi thẳng qua các lỗ tại vị trí Y = 0.5mm. Trên 
biểu đồ, quan sát thời điểm F̅ đạt giá trị tới hạn, 
tức khi mẫu bắt đầu bị phá hoại, cho thấy mẫu 1 
bắt đầu bị phá hoại sớm hơn mẫu 2. Chuyển vị 
cưỡng bức tương ứng với giá trị lớn nhất F̅ của 
mẫu 1 và 2 tương ứng là ua = 0.006mm, ub = 
0.0072mm. 
So sánh hai mẫu m3 và m4 (Hình 5), các 
lỗ được bố trí hai bên vết nứt ban đầu. Tương tự 
như ở mẫu m1 và m2, tọa độ theo phương X ở hai 
mẫu hoàn toàn giống nhau nhưng khác nhau theo 
phương Y. Khoảng cách giữa tâm các lỗ tương 
ứng là d = 0.08mm và d = 0.1386 mm. Quan sát 
phân bố của trường d và quan hệ lực - chuyển vị 
ở hai mẫu m3 và m4, các lỗ rỗng được bố trí hai bên vết nứt ban đầu, cho thấy: (i) phân bố của trường d 
là giống nhau, vết nứt tiến triển từ vết nứt ban đầu và đi thẳng cho đến khi mẫu bị phá hoại hoàn toàn; (ii) 
không có sự khác biệt giữa quan hệ lực chuyển vị (trên Hình 6 đường cong F̅ – u̅ của hai mẫu m3 và m4 
chồng khít lên nhau dù bố trí lỗ rỗng ở hai mẫu là khác nhau). Điều này cho thấy các lỗ rỗng ở vùng này 
không ảnh hưởng đến giá trị lực tới hạn cũng như thời điểm đạt giá trị lực này, đồng thời giá trị lớn nhất F̅ 
ở hai mẫu m3 và m4 cũng lớn hơn mẫu a và d.
5. Kết luận
Bài báo này dùng một mô hình trường pha được sử dụng để nghiên cứu sự ảnh hưởng phân bố của 
lỗ rỗng đến sự hình thành và phát triển của vết nứt. Kết quả bước đầu cho thấy: các lỗ rỗng bố trí ở 2 bên 
vết nứt ban đầu không ảnh hưởng đến sự khởi tạo và lan truyền vết nứt. Vết nứt không phụ thuộc vào sự 
phân bố của lỗ xung quanh vết nứt ban đầu mà chỉ phụ thuộc vào môi trường phía trước vết nứt, vùng mà 
vết nứt được dự kiến phát triển. Kết quả này phù hợp với sự phân bố ứng suất ở hai bên vết nứt ban đầu 
và phù hợp với kết luận rằng độ cứng phá hủy chịu ảnh hưởng lớn của sự phân bố lỗ rỗng xung quanh vết 
nứt [21]. Trong bài báo cũng đã xây dựng lại công thức phần tử hữu hạn để tìm trường hư hại của mô hình 
trường pha áp dụng cho bài toán phá hủy, là một mô hình mới, nhận được nhiều sự quan tâm trong lĩnh vực 
cơ học phá hủy hiện nay. Ảnh hưởng của khoảng cách, mật độ lỗ rỗng phía trước vết nứt tới sự phát triển 
vết nứt sẽ là chủ để nghiên cứu tiếp theo của nhóm tác giả.
Lời cảm ơn: Nghiên cứu này được tài trợ bởi Quỹ phát triển Khoa học và Công nghệ quốc gia (NAFOST-
ED) trong đề tài mã số 107.02-2017.309.
Tài liệu tham khảo
1. Ngo D., Scordelis A.C. (1967), “Finite element analysis of reinforced concrete beams”, Journal of Ameri-
can Concrete Institute, 64:152-163. 
2. Ingraffea A., Saouma V. (1984), Numerical modelling of discrete crack propagation in reinforced and plain 
concrete, In: Sih G, Di Tommaso A, editors, Fracture mechanics of concrete, Dordrecht: Martinus Nijhoff 
Publishers; 171-225.
3. Camacho G.T., Ortiz M. (1996), “Computational modeling of impact damage in brittle materials”, Interna-
tional Journal of Solids and Structures, 33:1267-1282.
4. Belytschko T., Black T. (1999), “Elastic crack growth in finite elements with minimal remeshing”, Interna-
tional Journal for Numerical Methods in Engineering, 45:601-20.
5. Daux C., Moës N., Dolbow J., Belytschko T. (2000), “Arbitrary branched and intersecting cracks with the 
extended finite element method”, International Journal for Numerical Methods in Engineering; 48:1741-60.
6. Moës N., Dolbow J., Belytschko T. (1999), “A finite element method for crack growth without remeshing”, 
International Journal for Numerical Methods in Engineering, 46(1):131-56.
Hình 6. So sánh các mối quan hệ lực (F) - chuyển vị (u) 
của các mẫu
107TẬP 11 SỐ 509 - 2017
KẾT QUẢ NGHIÊN CỨU VÀ ỨNG DỤNG
7. Moran B., Belytschko T., Sukumar N., Moës N. (2000), “Extended finite element method for three-dimen-
sional crack modeling”, International Journal for Numerical Methods in Engineering; 48(11):1549-70. 
8. Aranson I.S., Kalatsky V.A., Vinokur V.M. (2000), “Continuum field description of crack propagation”, 
Physical Review Letters, 85:118-121.
9. Karma A., Kessler D.A., Levine H. (2001), “Phase-field model of mode III dynamic fracture”, Physical 
Review Letters, 87:045501.
10. Hakim V., Karma A. (2009), “Laws of crack motion and phase-field models of fracture”, Journal of the 
Mechanics and Physics of Solids, 57:342-368.
11. Francfort G.A., Marigo J.J. (1998), “Revisiting brittle fractures as an energy minimization problem”, Journal 
of the Mechanics and Physics of Solids, 46:1319-1342.
12. Bourdin B., Francfort G.A., Marigo J.J. (2000), “Numerical experiments in revisited brittle fracture”, Journal 
of the Mechanics and Physics of Solids, 48:797-826.
13. Bourdin B., Francfort G. A., Marigo J.J. (2008), “The variational approach to fracture”, Journal of Elasticity, 
91:5-148.
14. Larsen C.J., Ortner C., Süli E. (2010), “Existence of solutions to a regularized model of dynamic fracture”, 
Mathematical Models and Methods in Applied Sciences, 20:1021-1048.
15. Bourdin B., Larsen C.J., Richardson C. (2011), “A time-discretemodel for dynamic fracture based on 
crack regularization”, International Journal of Fracture, 168:133-14.
16. Schlüter A., Willenbücher A., Kuhn C., Müller R. (2014), “Phase field approximation of dynamic brittle 
fracture”, Computational Mechanics, 54:1141-1161.
17. Verhoosel C.V., de Borst R. (2013), “A phase-field model for cohesive fracture”, International Journal for 
Numerical Methods in Engineering, 96:43-62.
18. May S., Vignollet J., Borst R.D.(2015), “A numerical assessment of phase models for brittle and cohesive 
fracture: G-convergence and stress oscilla-tions”, European Journal of Mechanics-A/Solids, 52:72-84.
19. Miehe C., Hofacker M., Welschinger F. (2010), “A phase model for rate-independentcrack propagation: 
Robust algorithmic implementation based on operator splits”, Computer Methods in Applied Mechanics and 
Engineering, 199(45):2765-2778.
20. Nguyen T.T., Yvonnet J., Zhu Q-Z., Bornert M., Chateau C. (2015), “A phase method to simulatecrack 
nucleation and propagation in strongly heteroge-neous materials from direct imaging of their microstructure”, 
Engineering Fracture Mechanics, 139:18-39.
21. Cramer M., Sevostianov M. (2009), “Effect of pore distribution on elastic stiffness and fracture toughness 
of porous media”, International Journal of Fracture, doi:10.1007/s10704-009-9416-6.