Help! I want to use LAMMPS to create a simulation for EDM, but I'm encountering some problems

Fykon · May 14, 2024, 12:06pm

Hi !
I am a first-year graduate student, recently began to use LAMMPS application in EDM analysis, the specific effect I want to do and my code is as follows, please help me to see what the problem is.

Initialization

units metal

dimension 3

boundary p p p

atom_style atomic

#Define the simulation area

region whole block 0 300 0 300 0 500 units box

create_box 2 whole

Create an electrode

lattice fcc 3.615

region electrode block 100 200 100 150 230 330 units box

create_atoms 1 region electrode

group electrode region electrode

Create artifacts

lattice fcc 3.615

region workpiece block 100 200 100 150 100 200 units box

create_atoms 2 region workpiece

group workpiece region workpiece

Define the side boundary layer area of the electrode

region el_side2 block 195 200 100 150 230 325 units box

region el_side3 block 105 195 100 105 230 325 units box

region el_side4 block 105 195 145 150 230 325 units box

region el_side5 block 100 200 100 150 325 330 units box

region wp_side1 block 100 105 100 150 100 200 units box

region wp_side2 block 195 200 100 150 100 200 units box

region wp_side3 block 105 195 100 105 100 200 units box

region wp_side4 block 105 195 145 150 100 200 units box

group el_side2 region el_side2

group el_side3 region el_side3

group el_side4 region el_side4

group el_side5 region el_side5

group wp_side1 region wp_side1

group wp_side2 region wp_side2

group wp_side3 region wp_side3

group wp_side4 region wp_side4

The thermostatic region of the electrode

region el_hwc block 105 195 105 145 305 325 units box

region wp_hwc block 105 195 105 145 100 120 units box

region hwc union 2 el_hwc wp_hwc

group el_hwc region el_hwc

group wp_hwc region wp_hwc

group hwc region hwc

Newtonian layer

region el_ndc block 105 195 105 145 230 305 units box

region wp_ndc block 105 195 105 145 180 200 units box

region ndc union 2 el_ndc wp_ndc

group wp_ndc region wp_ndc

group el_ndc region el_ndc

group ndc region ndc

Define atomic type and mass

mass 1 63.546 # Cu

mass 2 63.546 # Cu

Select the applicable force field

pair_style eam/alloy

pair_coeff * * Cu_Zhou04.eam.alloy Cu Cu

Fixed the motion of atoms in the boundary layer on four sides

fix fix_el_side2 el_side2 setforce 0.0 0.0 0.0

fix fix_el_side4 el_side4 setforce 0.0 0.0 0.0

fix fix_el_side5 el_side5 setforce 0.0 0.0 0.0

fix fix_wp_side1 wp_side1 setforce 0.0 0.0 0.0

fix fix_wp_side2 wp_side2 setforce 0.0 0.0 0.0

fix fix_wp_side4 wp_side4 setforce 0.0 0.0 0.0

#Set the initial temperature

velocity all create 300 12345

Application of constant temperature layer

fix nvt_hwc hwc nvt temp 300 300 0.1

Workpiece temperature

fix 2 wp_ndc nvt temp 300 300 0.1

Apply temperature control to the electrode

fix 3 el_ndc nvt temp 300 300 0.1

A variable that controls whether the heat source is activated or not, stop heating after 15000 steps

variable heat_active equal “step < 15000”

variable heat_mult equal “v_heat_active * 20.8”

variable heat_mult_2 equal “v_heat_active * 8.8”

The first heat source

variable x0 equal 120

variable y0 equal 115

variable z0 equal 240 # Central position of heat source

variable Pm equal ${heat_mult} # Maximum heating intensity

variable k equal 0.03 # Parameters that control the range of Gaussian distribution

variable fx atom v_Pmexp(-v_k((x-v_x0)^2+(y-v_y0)^2+(z-v_z0)^2))

variable fy atom v_Pmexp(-v_k((x-v_x0)^2+(y-v_y0)^2+(z-v_z0)^2))

variable fz atom v_Pmexp(-v_k((x-v_x0)^2+(y-v_y0)^2+(z-v_z0)^2))

Applying the force of Gaussian distribution to the Newtonian layer

fix heat el_ndc addforce v_fx v_fy v_fz

The second heat source

variable x0_2 equal 155

variable y0_2 equal 130

variable z0_2 equal 210

variable Pm2 equal ${heat_mult_2}

variable k equal 0.03

variable fx_2 atom v_Pm2exp(-v_k((x-v_x0_2)^2+(y-v_y0_2)^2+(z-v_z0_2)^2))

variable fy_2 atom v_Pm2exp(-v_k((x-v_x0_2)^2+(y-v_y0_2)^2+(z-v_z0_2)^2))

variable fz_2 atom v_Pm2exp(-v_k((x-v_x0_2)^2+(y-v_y0_2)^2+(z-v_z0_2)^2))

fix heat2 wp_ndc addforce v_fx_2 v_fy_2 v_fz_2

thermo 500

thermo_style custom step temp pe ke etotal press pxx pyy pzz

dump 1 all custom 500 dump.Cu.xyz id type x y z

run 60000

I follow the production of this article, I feel that there is something wrong with my Gaussian heat source setting, because this article does not specify how to set it.(Influence of discharge gap on material removal and melt pool movement in EDM discharge process    (https://doi.org/10.1007/s00170-021-08577-z))

I think the key lies in the setting of the heat source. I am not sure how to set the heat source properly to achieve the effect in the article. Of course, my other settings are also rough. Please give suggestions.

hothello · May 14, 2024, 9:18pm

The first problem is that you didn’t read the forum guidelines, but you can still do it and edit your first post accordingly.

The second problem is that you need to communicate your problem effectively, describing rigorously what problem(s) you encountered, providing context, etc. It is way too vague to state that “the simulation doesn’t work”, as no one is interested in doing the homework on your behalf.

Regarding your script, I have made the following changes:

Simplified the definition of fixed and mobile regions,
Used the atom types to distinguish between fixed and mobile atoms,

image825×768 98.8 KB
Changed the position of the heat source (please check z2),
Added a random component to fix addforce to mimic the effect of a heat source and changed the Gaussian definition to my own taste,
Removed the fix setforce, as it suffices to integrate the mobile atoms.
Thermalised the mobile atoms before applying the heat source,
Used an NVE integrator for the atoms on which the force is applied,

Here is the revised script, with proper quotation:

# Initialization
units           metal
dimension       3
boundary        p p p
atom_style      atomic
processors      3 2 2

# Define the simulation area
region whole block 0 300 0 300 0 500 units box
create_box 2 whole

# Create an electrode
lattice fcc 3.615
region electrode block 100 200 100 150 230 330 units box
create_atoms 1 region electrode
group electrode region electrode

# Create artifacts
lattice fcc 3.615
region workpiece block 100 200 100 150 100 200 units box
create_atoms 1 region workpiece
group workpiece region workpiece

# Define atomic type and mass
mass            1    63.546 # Cu
mass            2    63.546 # Cu

# Select the applicable force field
pair_style eam/alloy
pair_coeff * * Cu_zhou.eam.alloy Cu Cu

# Define the side boundary layer area of the electrode
region          el_thermo block 105 195 105 145 230 325 units box
region          el_side   intersect 2 electrode el_thermo
region          wp_thermo block 105 195 105 145 105 200 units box
region          wp_side   intersect 2 workpiece wp_thermo
region          mobile    union 2 el_thermo wp_thermo
region          side      union 2 el_side wp_side
group           mobile    region mobile
group           side      region side
set group side type 2

# The thermostating region of the electrode
region          el_hwc block 105 195 105 145 250 325 units box
region          wp_hwc block 105 195 105 145 105 180 units box
region          hwc    union 2 el_hwc wp_hwc # nvt
group           hwc    region hwc

# Newtonian layer
region          el_ndc intersect 2 el_thermo el_hwc # addforce
region          wp_ndc intersect 2 wp_thermo wp_hwc # addforce
region          ndc union 2 el_ndc wp_ndc # nve
group           wp_ndc region wp_ndc
group           el_ndc region el_ndc
group           ndc region ndc

# Thermodynamic output.
comm_style      tiled
balance         1.0 rcb
compute         mobile_temp mobile temp
thermo          3000
thermo_style    custom step temp pe ke etotal press pxx pyy pzz
thermo_modify   temp mobile_temp
dump            1        all      custom 500 dump.Cu.xyz id type x y z

# First, relax the Cu blocks.
velocity        all create 300 68748
fix             NVT      mobile      nvt temp 300 300 0.1
run             10000
unfix           NVT

# A variable that controls whether the heat source is activated or not, stop heating after 15000 steps
variable heat_active equal "step < 15000"
variable heat_mult equal "v_heat_active * 20.8"
variable heat_mult_2 equal "v_heat_active * 8.8"

# The first heat source
variable x0 equal 120
variable y0 equal 115
variable z0 equal 240 # Central position of heat source
variable Pm equal ${heat_mult} # Maximum heating intensity
variable k equal 0.03 # Parameters that control the range of Gaussian distribution
variable fx atom v_Pm*exp(-v_k*(x-v_x0)^2)*random(-1,1,37324)
variable fy atom v_Pm*exp(-v_k*(y-v_y0)^2)*random(-1,1,82396)
variable fz atom v_Pm*exp(-v_k*(z-v_z0)^2)*random(-1,1,72383)

# Applying the force of Gaussian distribution to the Newtonian layer
fix             heat     el_ndc   addforce v_fx v_fy v_fz every 10

# The second heat source
variable x2 equal 155
variable y2 equal 130
variable z2 equal 190
variable Pm2 equal ${heat_mult_2}
variable fx_2 atom v_Pm2*exp(-v_k*(x-v_x2)^2)*random(-1,1,24535)
variable fy_2 atom v_Pm2*exp(-v_k*(y-v_y2)^2)*random(-1,1,61699)
variable fz_2 atom v_Pm2*exp(-v_k*(z-v_z2)^2)*random(-1,1,12293)
fix             heat2    wp_ndc   addforce v_fx_2 v_fy_2 v_fz_2 every 10

# Integration and temperature control on the electrodes' bulk.
fix             NVT      hwc      nvt temp 300 300 0.1

# Integration only on the electrodes' surface.
fix             NVE      ndc      nve
run             60000

WARNING: One or more atoms are time-integrated more than once. I am not sure why this happens, as the regions hwc and ndc should be defined as not to have any overlap. I leave the debugging to you (or any other curious user).

hothello · May 16, 2024, 9:19am

Please check this discussion on region intersect. I am attaching the revised script, which runs without errors. Your homework is to define a fix to compute the amount of energy added to the two Cu blocks over time: with your initial choice (and the NVE integrator for the heated region), the surfaces kind of evaporate at supersonic speed. Instead, you may want to control carefully the heat flux added to the system.
The new sample (cut into half to visualise the internal structure):

# Initialization
units           metal
dimension       3
boundary        p p p
atom_style      atomic
#processors      3 2 2

# Define the simulation area
region whole block 0 300 0 300 0 500 units box
create_box 2 whole

# Create two electrodes
lattice         fcc 3.615
region          electrode block 100 200 100 150 230 330 units box
region          workpiece block 100 200 100 150 100 200 units box
create_atoms    1 region electrode
create_atoms    1 region workpiece
group           electrode region electrode
group           workpiece region workpiece

# Define atomic type and mass
mass            1    63.546 # Cu
mass            2    63.546 # Cu

# Select the applicable force field
pair_style eam/alloy
pair_coeff * * Cu_zhou.eam.alloy Cu Cu

# Partition the electrodes into fixed and mobile regions.
region          el_bulk   block 105 195 104 145 230 327 units box
region          el_surf   block 100 200 100 150 230 235 units box
region          el_thermo union 2 el_bulk el_surf
region          wp_bulk   block 105 195 104 145 104 200 units box
region          wp_surf   block 100 200 100 150 195 200 units box
region          wp_thermo union 2 wp_bulk wp_surf
group           el_thermo region el_thermo
group           wp_thermo region wp_thermo
group           el_side subtract electrode el_thermo
group           wp_side subtract workpiece wp_thermo
group           mobile union el_thermo wp_thermo
group           side union el_side wp_side
set             group side type 2

# The thermostating region of the electrode
region          el_hwc block 105 195 104 145 250 327 units box
region          wp_hwc block 105 195 104 145 104 180 units box
group           el_hwc region el_hwc
group           wp_hwc region wp_hwc
group           hwc    union el_hwc wp_hwc # nvt

# Newtonian layer
group           el_ndc subtract el_thermo el_hwc # addforce
group           wp_ndc subtract wp_thermo wp_hwc # addforce
group           ndc union el_ndc wp_ndc # nve

# Thermodynamic output.
comm_style      tiled
balance         1.0 rcb
compute         mobile_temp mobile temp
thermo          500
thermo_style    custom step temp pe ke etotal press pxx pyy pzz
thermo_modify   temp mobile_temp flush yes
dump            1        all      custom 500 dump.Cu.xyz id type x y z

# First, relax the Cu blocks.
velocity        all create 300 68748
fix             NVT      mobile      nvt temp 300 300 0.1
run             3000
unfix           NVT

# A variable that controls whether the heat source is activated or not, stop heating after 15000 steps
variable heat_active equal "step < 15000"
variable heat_mult  equal "v_heat_active * 10."
variable heat_mult2 equal "v_heat_active * 5."

# The first heat source
variable        x1 equal 120
variable        y1 equal 115
variable        z1 equal 240 # Central position of heat source
variable        Pm equal ${heat_mult} # Maximum heating intensity
variable        k equal 0.03 # Parameters that control the range of Gaussian distribution
variable fx atom v_Pm*exp(-v_k*(x-v_x1)^2)*random(-1,1,37324)
variable fy atom v_Pm*exp(-v_k*(y-v_y1)^2)*random(-1,1,82396)
variable fz atom v_Pm*exp(-v_k*(z-v_z1)^2)*random(-1,1,72383)

# Applying the force of Gaussian distribution to the Newtonian layer
fix             heat     el_ndc   addforce v_fx v_fy v_fz every 100

# The second heat source
variable        x2 equal 155
variable        y2 equal 130
variable        z2 equal 190
variable Pm2 equal ${heat_mult2}
variable fx2 atom v_Pm2*exp(-v_k*(x-v_x2)^2)*random(-1,1,24535)
variable fy2 atom v_Pm2*exp(-v_k*(y-v_y2)^2)*random(-1,1,61699)
variable fz2 atom v_Pm2*exp(-v_k*(z-v_z2)^2)*random(-1,1,12293)
fix             heat2    wp_ndc   addforce v_fx2 v_fy2 v_fz2 every 100

# Integration and temperature control on the electrodes' bulk.
fix             NVT      hwc      nvt temp 300 300 0.1

# Integration only on the electrodes' surface.
fix             NVE      ndc      nve
run             60000

Fykon · May 18, 2024, 4:58am

Thank you for your valuable advice and assistance. I have just come across your response, and I will carefully analyze it.

Fykon · May 18, 2024, 4:59am

After that I tried to use the speed of giving the atoms a Gaussian heat source structure, and it seemed to get close to the result I wanted