LAMMPS Fracture

From EVOCD
(Difference between revisions)
Jump to: navigation, search
(Abstract)
(LAMMPS input script)
Line 24: Line 24:
 
|<pre>
 
|<pre>
 
############################################################################
 
############################################################################
# Interfacial fracture
+
# Interfacial Fracture
# Mark Tschopp, 2010
+
# Mark Tschopp, Nathan Rhodes 2011
  
# lmp_exe -var datfile Fe_110_sig3.txt -var strain 0.001 -var nloop 100 -var minlength 20 < in.gb_fracture.txt
+
# lmp_exe -var datfile test.data -var nloop 100 < in.gb_fracture.txt
 +
 
 +
# Simulation deletes atoms outside of +/- deldist from GB and constrains and pulls
 +
# atoms outside of +/- fixdist from GB to fracture the GB
 
############################################################################
 
############################################################################
  
variable datfile index Fe_110_sig3.txt
+
#variable datfile index test.data
 
variable strain equal 0.001
 
variable strain equal 0.001
variable nloop equal 100
+
#variable nloop equal 100
#variable repl equal 1
+
variable repx equal 1
 +
variable repz equal 1
 
variable strain2 equal "1+v_strain"
 
variable strain2 equal "1+v_strain"
 +
variable deldist equal 50
 +
variable fixdist equal 45
  
 
######################################
 
######################################
Line 58: Line 64:
 
print "lz: ${zlen}"
 
print "lz: ${zlen}"
  
# Determine number of increments for displacement grid in the in-plane GB directions  
+
# Determine number of increments for displacement grid in the in-plane GB directions
variable xrepl equal "ceil(v_minlength / v_xlen)"  
+
######### This is outdated, as replicated directions are now manually adjusted ##########
variable zrepl equal "ceil(v_minlength / v_zlen)"  
+
#variable xrepl equal "ceil(v_minlength / v_xlen)"  
 +
#variable zrepl equal "ceil(v_minlength / v_zlen)"  
  
replicate ${xrepl} 1 ${zrepl}
+
#OLD================================================================replicate ${xrepl} 1 ${zrepl}
 +
replicate ${repx} 1 ${repz}
  
 
######################################
 
######################################
Line 69: Line 77:
 
pair_coeff * * Fe-C_Hepburn_Ackland.eam.fs Fe C
 
pair_coeff * * Fe-C_Hepburn_Ackland.eam.fs Fe C
  
 +
 +
compute stress all stress/atom
 +
compute stress1 all reduce sum c_stress[1]
 +
compute stress2 all reduce sum c_stress[2]
 +
compute stress3 all reduce sum c_stress[3]
 +
compute stress4 all reduce sum c_stress[4]
 +
compute stress5 all reduce sum c_stress[5]
 +
compute stress6 all reduce sum c_stress[6]
 
##########################################
 
##########################################
 
# Minimize first?
 
# Minimize first?
 
reset_timestep 0
 
reset_timestep 0
 
thermo 10
 
thermo 10
thermo_style custom step lx ly lz press pxx pyy pzz pe
+
thermo_style custom step lx ly lz press pxx pyy pzz pe c_stress1 c_stress2 c_stress3 c_stress4 c_stress5 c_stress6
 
min_style cg
 
min_style cg
 
fix 1 all box/relax x 0.0 z 0.0 couple none vmax 0.001  
 
fix 1 all box/relax x 0.0 z 0.0 couple none vmax 0.001  
Line 93: Line 109:
 
variable p8 equal "pe"
 
variable p8 equal "pe"
  
fix equil1 all print 1 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8}" file data.${datfile}_${minlength}.txt screen no
+
#OLD================================================================fix equil1 all print 1 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8}" file data.${datfile}_${minlength}.txt screen no
 +
fix equil1 all print 1 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8}" file data.${datfile}.txt screen no
 
fix 1 all nve
 
fix 1 all nve
 
run 1
 
run 1
Line 100: Line 117:
 
variable pressf equal ${pressf1}
 
variable pressf equal ${pressf1}
  
 +
##########################################
 +
# Create cfg files with stress in y direction for AtomEye viewing
 +
 +
dump 1 all cfg 500 ${datfile}_*.cfg id type xs ys zs c_stress[2] fx fy fz
 +
#dump_modify 1 element Fe
 +
 +
 +
##########################################
 +
# Create regions for fixed boundary condition (FOR NEW METHOD)
 +
 +
region rlow block 0 200 -200 -${deldist} 0 200 units box
 +
region rhigh block 0 200 ${deldist} 200 0 200 units box
 +
region rgblow block 0 200 -50 -${fixdist} 0 200 units box
 +
region rgbhigh block 0 200 ${fixdist} 50 0 200 units box
 +
 +
group glow region rlow
 +
group ghigh region rhigh
 +
group gbhigh region rgbhigh
 +
group gblow region rgblow
 +
 +
# Delete atoms far from grain boundary
 +
 +
delete_atoms group glow
 +
delete_atoms group ghigh
 +
 +
# Put fixed boundary condition on edge atoms by setting forces to zero
 +
fix 2 gbhigh setforce 0 0 0
 +
fix 3 gblow setforce 0 0 0
 
##########################################
 
##########################################
 
# MS Deformation loop
 
# MS Deformation loop
Line 107: Line 152:
  
 
# Increase box bound and minimize again
 
# Increase box bound and minimize again
reset_timestep 0
+
#reset_timestep 0
 
#displace_box all y scale ${strain2}
 
#displace_box all y scale ${strain2}
 
#fix 1 all box/relax x 10000.0 z 10000.0 couple none vmax 0.001  
 
#fix 1 all box/relax x 10000.0 z 10000.0 couple none vmax 0.001  
displace_box all y delta -${lydelta} ${lydelta} units box
+
 
 +
###########ORIGINAL METHOD#############
 +
#displace_box all y delta -${lydelta} ${lydelta} units box
 +
 
 +
##############NEW METHOD###############
 +
displace_box gblow y delta -${lydelta} 0 units box
 +
displace_box gbhigh y delta 0 ${lydelta} units box
 +
 
 +
#displace_atoms gblow move 0 -${lydelta} 0 units box
 +
#displace_atoms gbhigh move 0 ${lydelta} 0 units box
 +
 
 
minimize 1.0e-25 1.0e-25 1000 10000
 
minimize 1.0e-25 1.0e-25 1000 10000
  
Line 119: Line 174:
 
print "Pressf: ${pressf}"
 
print "Pressf: ${pressf}"
 
print "Pdiff: ${pdiff}"
 
print "Pdiff: ${pdiff}"
if ${pdiff} > 10000 then "exit"
+
#if ${pdiff} > 10000 then "exit"
 
variable pressf1 equal pyy
 
variable pressf1 equal pyy
 
variable pressf equal ${pressf1}
 
variable pressf equal ${pressf1}

Revision as of 14:05, 10 May 2011

Contents

Abstract

This example shows how to run an atomistic simulation of fracture of an iron symmetric tilt grain boundary. A parallel molecular dynamics code, LAMMPS[1], is used to calculate stresses at the grain boundary as the strain of the bicrystal is incrementally increased. Matlab is used to plot a stress-strain curve, and Atomeye is used to visualize the simulation.

Author(s): Mark A. Tschopp, Nathan R. Rhodes

Corresponding Author: Mark Tschopp

Input

Description of Simulation

This molecular dynamics simulation calculates the stress-strain relationship of an iron symmetric tilt grain boundary under fracture. The grain boundary structure used in this example is a <100> Σ5(210) symmetric tilt grain boundary. The potential used to generate the structure, the Hepburn and Ackland (2008) Fe-C interatomic potential[2], is also used in this script. The simulation cell is defined such that the bicrystal is pulled in the y-direction, or perpendicular to the boundary interface, to increase strain. The strain in increased for a specified number of times in a loop, and the stress is calculated at each point before the simulation loops. The stress and strain values are output to a separate file which can be imported in a graphing application for plotting.

Grain boundary structure file

The grain boundary structure that was generated prior to this example can be found here. Store the text in "Fe_110_sig3.txt" to use it.

LAMMPS input script

This input script was run using the November 2010 version of LAMMPS. Changes in some commands in more recent versions may require revision of the input script. This script runs the simulation with a previously generated grain boundary file, which is fed to the variable "datfile." The variable "nloop" defines how many times the strain will be increased and the number of points at which the stress is calculated. To run this script, store it in "in.gb_fracture.txt" and use "lmp_exe < in.gb_fracture.txt" in a UNIX environment where "lmp_exe" refers to the LAMMPS executable.

############################################################################
# Interfacial Fracture
# Mark Tschopp, Nathan Rhodes 2011

# lmp_exe -var datfile test.data -var nloop 100 < in.gb_fracture.txt

# Simulation deletes atoms outside of +/- deldist from GB and constrains and pulls
# atoms outside of +/- fixdist from GB to fracture the GB
############################################################################

#variable datfile index test.data
variable strain equal 0.001
#variable nloop equal 100
variable repx equal 1
variable repz equal 1
variable strain2 equal "1+v_strain"
variable deldist equal 50
variable fixdist equal 45

######################################
# INITIALIZATION
units 		metal
dimension		3
boundary		p	p	p
atom_style		atomic
atom_modify map array

######################################
# SIMULATION CELL VARIABLES (in Angstroms)

read_data ${datfile}

#variable minlength equal 100
variable xlen equal lx
variable ylen equal ly
variable zlen equal lz

print "lx: ${xlen}"
print "ly: ${ylen}"
print "lz: ${zlen}"

# Determine number of increments for displacement grid in the in-plane GB directions
######### This is outdated, as replicated directions are now manually adjusted ##########
#variable xrepl equal "ceil(v_minlength / v_xlen)" 
#variable zrepl equal "ceil(v_minlength / v_zlen)" 

#OLD================================================================replicate ${xrepl} 1 ${zrepl}
replicate ${repx} 1 ${repz}

######################################
# INTERATOMIC POTENTIAL
pair_style	eam/fs
pair_coeff	* * Fe-C_Hepburn_Ackland.eam.fs Fe C


compute stress all stress/atom
compute stress1 all reduce sum c_stress[1]
compute stress2 all reduce sum c_stress[2]
compute stress3 all reduce sum c_stress[3]
compute stress4 all reduce sum c_stress[4]
compute stress5 all reduce sum c_stress[5]
compute stress6 all reduce sum c_stress[6]
##########################################
# Minimize first?
reset_timestep 0
thermo		10
thermo_style custom step lx ly lz press pxx pyy pzz pe c_stress1 c_stress2 c_stress3 c_stress4 c_stress5 c_stress6
min_style cg
fix 1 all box/relax x 0.0 z 0.0 couple none vmax 0.001 
minimize 1.0e-25 1.0e-25 1000 10000
unfix 1

variable ly1 equal ly
variable ly0 equal ${ly1}
variable lydelta equal "v_strain*v_ly0/2"

# Setup file output (time in ps, pressure in GPa)
variable p1 equal "(ly-v_ly0)/v_ly0"
variable p2 equal "-pxx/10000"
variable p3 equal "-pyy/10000"
variable p4 equal "-pzz/10000"
variable p5 equal "-pxy/10000"
variable p6 equal "-pxz/10000"
variable p7 equal "-pyz/10000"
variable p8 equal "pe"

#OLD================================================================fix equil1 all print 1 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8}" file data.${datfile}_${minlength}.txt screen no
fix equil1 all print 1 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8}" file data.${datfile}.txt screen no
fix 1 all nve
run 1
unfix 1
variable pressf1 equal pyy
variable pressf equal ${pressf1}

##########################################
# Create cfg files with stress in y direction for AtomEye viewing

dump 1 all cfg 500 ${datfile}_*.cfg id type xs ys zs c_stress[2] fx fy fz
#dump_modify 1 element Fe


##########################################
# Create regions for fixed boundary condition (FOR NEW METHOD)

region rlow block 0 200 -200 -${deldist} 0 200 units box
region rhigh block 0 200 ${deldist} 200 0 200 units box
region rgblow block 0 200 -50 -${fixdist} 0 200 units box
region rgbhigh block 0 200 ${fixdist} 50 0 200 units box

group glow region rlow
group ghigh region rhigh
group gbhigh region rgbhigh
group gblow region rgblow

# Delete atoms far from grain boundary

delete_atoms group glow
delete_atoms group ghigh

# Put fixed boundary condition on edge atoms by setting forces to zero
fix 2 gbhigh setforce 0 0 0
fix 3 gblow setforce 0 0 0
##########################################
# MS Deformation loop

variable a loop ${nloop}
label loop

# Increase box bound and minimize again
#reset_timestep 0
#displace_box all y scale ${strain2}
#fix 1 all box/relax x 10000.0 z 10000.0 couple none vmax 0.001 

###########ORIGINAL METHOD#############
#displace_box all y delta -${lydelta} ${lydelta} units box

##############NEW METHOD###############
displace_box gblow y delta -${lydelta} 0 units box
displace_box gbhigh y delta 0 ${lydelta} units box

#displace_atoms gblow move 0 -${lydelta} 0 units box
#displace_atoms gbhigh move 0 ${lydelta} 0 units box

minimize 1.0e-25 1.0e-25 1000 10000

run 1

print "Pressf: ${pressf}"
variable pdiff equal "pyy - v_pressf"
print "Pressf: ${pressf}"
print "Pdiff: ${pdiff}"
#if ${pdiff} > 10000 then "exit"
variable pressf1 equal pyy
variable pressf equal ${pressf1}

next a
jump in.gb_fracture.txt loop 

unfix equil1

######################################
# SIMULATION DONE
print "All done"

Output

LAMMPS datafile

The following file, named "data.Fe_100_sig52_10.txt" should have been created in addition to the log.lammps file. This file stores strain information in the first column, stress tensor information in the second through seventh columns, and stores the total potential energy of the cell in the eight column. The simulation should have looped 100 times (as per the "nloop" variable), so there should be 100 entries (which end at a strain of 0.1) plus the initial entry of stress and strain at zero.

# Fix print output for fix equil1
0 -2.219884235e-05 -0.09436668009 5.407552011e-06 1.221556825e-05 1.110009611e-09 4.462377068e-12 -48716.72405
0.001 0.1468487572 0.1387359115 0.1431841846 -6.180770704e-09 -2.564170766e-12 -2.40791993e-12 -48716.76403
0.002 0.2598702618 0.4058469325 0.2847557426 -4.808171999e-12 9.566142713e-16 -3.41722654e-15 -48716.52319
0.003 0.3720478051 0.6731956939 0.4254156925 1.311052574e-10 3.449206559e-15 -2.636944941e-14 -48716.04581
0.004 0.483977058 0.9402201271 0.5652029445 -6.629325468e-12 1.129857588e-15 2.416966446e-14 -48715.3319
0.005 0.5948099239 1.207643528 0.7041151176 -9.516316287e-13 2.79361336e-15 7.471854238e-15 -48714.38144
...


Post-Processing

Stress-Strain Plot

The stress-strain curve in Figure 1 can be generated using the following MATLAB script. The "exportfig" command saves the plot to a tiff file, but the plot can also be saved as a Mathcad figure once it appears. The three principle stresses are plotted: pxx in red, pyy in blue, and pzz in green. The shear stresses (not plotted) are stored in the fourth, fifth, and sixth columns of the stress variable.

%% Plot stress-strain curve from single GB fracture datafile

clear all

d = dir('data.Fe_100STGB*_min.data_20.txt');
for i = 1%:length(d)
    fname = d(i).name;
    A = importdata(fname);
    %Define stress and strain
    strain = A.data(:,1);
    stress = A.data(:,2:7);
    %Plot data
    plot(strain,stress(:,1),'-dr','LineWidth',2,'MarkerEdgeColor','r',...
        'MarkerFaceColor','r','MarkerSize',10), hold on
    plot(strain,stress(:,2),'-ob','LineWidth',2,'MarkerEdgeColor','b',...
        'MarkerFaceColor','b','MarkerSize',5), hold on
    plot(strain,stress(:,3),'-og','LineWidth',2,'MarkerEdgeColor','g',...
        'MarkerFaceColor','g','MarkerSize',5), hold on
    axis square
    %ylim([0 10])
    set(gca,'LineWidth',2,'FontSize',24,'FontWeight','normal','FontName','Times')
    set(get(gca,'XLabel'),'String','Strain','FontSize',32,'FontWeight','bold',...
        'FontName','Times')
    set(get(gca,'YLabel'),'String','Stress (GPa)','FontSize',32,'FontWeight','bold',...
        'FontName','Times')
    set(gcf,'Position',[1 1 round(1000) round(1000)])

    %Export the figure to a tif file
    exportfig(gcf,strrep(fname,'_min.data','.tif'),'Format','tiff',...
        'Color','rgb','Resolution',300)
    close(1)
end
Figure 1. Stress-strain curves for fracture of a <100> Σ5(210) Fe symmetric tilt grain boundary. The x, y, and z principle stresses are plotted in red, blue, and green, respectively.

Go Back

Acknowledgments

The authors would like to acknowledge funding for this work through the Department of Energy.

References

  1. S. Plimpton, "Fast Parallel Algorithms for Short-Range Molecular Dynamics," J. Comp. Phys., 117, 1-19 (1995).
  2. D.J. Hepburn and G.J. Ackland, "Metallic-covalent interatomic potential for carbon in iron," Phys. Rev. B 78, 165115 (2008).
Personal tools
Namespaces

Variants
Actions
home
Materials
Material Models
Design
Resources
Projects
Education
Toolbox