Problems with fixed areas during tensile simulation

​Hello,

I am running a tensile simulation and facing some problems with my results. For the tensile simulation I defined some fixed areas and increased the box bound with a predefined strain via the “change_box” command.

Looking at the stress-strain curve, my result seems to be quit nice, but after the separation, there is still some stress in the fixed areas. So after the separation of the material I have an increasing strain, because of the stress left in the fixed areas.

I am not sure if this is a numerical lammps problem or I did something wrong. I tried a few ways, but nothing really worked out. Do you maybe have some suggestions for me to solve that problem?

Here is my Lammps Input File:

------------------------ INITIALIZATION ----------------------------

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

----------------------- ATOM DEFINITION ----------------------------

lattice fcc ${latparam}
region whole block 0 30 -30 30 0 5 units box
create_box 2 whole
region middle block INF INF INF INF INF INF units box
create_atoms 2 region middle
group middle type 1

------------------------ FORCE FIELDS ------------------------------

pair_style eam/alloy
pair_coeff * * Al99.eam.alloy Al Al

---------------------- COMPUTE SETTINGS ----------------------------

compute csym all centro/atom fcc
compute eng all pe/atom #stores in eng potenergy of each atom
compute eatoms all reduce sum c_eng
compute peratom all pe/atom
compute stress all stress/atom NULL
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 SYSTEM ----------------------------

reset_timestep 0
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_stress1 c_stress2 c_stress3 c_stress4 c_stress5 c_stress6 c_eatoms temp fmax
fix 1 all box/relax y 1000.0 vmax 0.01
min_style cg
minimize 1.0e-15 1.0e-15 100000 100000
unfix 1

------------------------ STRESS STRAIN OUTPUT ----------------------

variable strain equal 0.001
variable ly1 equal ly
variable ly0 equal {ly1} variable lydelta equal "v_strain*v_ly0/2" 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" fix equil1 all print 1 "{p1} {p2} {p3} {p4} {p5} {p6} {p7} ${p8}" file data.RESULT.txt screen no
fix 1 all nve
run 1
unfix 1

----------------------- TENSILE SIMULATION -------------------------

region rgblow block 0 200 -200 -29 0 200 units box
region rgbhigh block 0 200 29 200 0 200 units box
group gbhigh region rgbhigh
group gblow region rgblow
fix 2 gbhigh setforce 0 0 0
fix 3 gblow setforce 0 0 0
variable nloop equal 600
variable a loop {nloop} label loop change_box gblow y delta -{lydelta} 0 remap units box
change_box gbhigh y delta 0 ${lydelta} remap units box
minimize 1.0e-35 1.0e-35 10000 100000
run 1
next a
jump single_tensile.in loop
unfix equil1

You can have a look at my results following this link:

Best Regards
Philipp

What is a fixed area? Why would you expect to

deform the box, and then have the stress everywhere

be 0.0? A more typical way to map a stress/strain curve

is to use fix deform during a running simulation.

Steve

Hi Steve,

thank you for your reply. I defined a fixed area where I set the forces down to zero. I expect that the stress will be zero after the separation of the material, as it is in my non fixed area. I tried the fix deform command, but it didn´t worked out. For tensile simulation in y-direction I expected that the max stress value is independent from the model size in y- direction. So I ran a simulation for 4 different models with different sizes in y-direction and received different values for the max stress value, using the fix deform command. Be defining a fixed area and working with the change_box command a got my expected result. Running a simulation for 4 different modelsizes in y-direction a received the same max stress values. The only problem I was facing, that the stress increases after the separation of the material.

Philipp

ok - if you’re satisfied, great. I don’t understand the connection

between setting force to 0 and expecting the stress to be 0.

In LAMMPS, the virial (stress) is calculated in the pair interactions.

There is no way to zero it. The fix setforce command operates

after all the pairwise interaction and simply zeroes the forces. It

will not zero the virial.

Steve

If I understand your script correctly, you have a 3D periodic cell containing a perfect crystal at zero temperature that you are incrementally stretching in the y direction using repeated change_box commands with remap. After each increment you do a minimize to relax the atom positions while freezing the two layers at ymax and ymin.

There are a lot of issues with this set up, mainly related to the fact that you have eliminated temperature, surfaces, defects, so you have suppressed all the normal tensile failure nucleation mechanisms. I think the crystal will preserve perfect symmetry up to the separation point and will only separate when the potential energy crosses a maximum w.r.t. Ly, which probably will result in a yield strength that is unphysically large. Once that happens, the minimization will generate a somewhat arbitrary cutting plane and all the atoms will relax very abruptly into two layers on opposites sides of the two layers of fixed atoms at Ymin and Ymax. Note that because you have periodicity in y, these two layers will behave like a single layer that just happens to have some frozen atoms in the middle. Probably the ever increasing lattice spacing in y for the two layers of frozen atoms is what is causing your tensile stress to keep increasing.

That said, if you want to see the stress go to zero, do the following:

  1. Turn off periodic boundaries in the y direction, so you have an infinite slab geometry
  2. After the slab separates into two slabs, exit the loop
  3. unfix the setforce fixes
  4. Do another minimization

At this point, since all the forces are zero, the stress will be zero.

As Steve suggested, a better way to run this simulation would be with fix deform. When you used it before, you were probably getting a yield stress that depended on sample length because:

  1. The frozen atoms in the layers at Ymin and Ymax were being included in the reported stress
  2. You were getting interactions between atoms across the periodic boundary
  3. There may have been other problems in your setup too

In order to clearly eliminate 1 and 2, I suggested turning off periodicity in y, and obtaining the average stress in the sample by averaging the per-atom stress of atoms that are more than the force cutoff from the boundary layers. You should also consider:

  1. Adding temperature
  2. Adding free surfaces normal to x or z directions
  3. Adding a notch in the free surfaces
  4. Reading the prior literature on MD simulation of tensile failure

Thank you for your reply. It really helped me and I am thankful about your suggestions. I am working on it.
Best Regards and again Thank you all!

Philipp