forked from lijiext/lammps
140 lines
6.5 KiB
Plaintext
140 lines
6.5 KiB
Plaintext
####################################################################################################
|
|
#
|
|
# TLSPH example: elongate a 2d strip of aluminum py pulling its ends apart
|
|
#
|
|
# unit sytem: GPa / mm / ms
|
|
#
|
|
####################################################################################################
|
|
|
|
|
|
####################################################################################################
|
|
# MATERIAL PARAMETERS
|
|
####################################################################################################
|
|
variable E equal 70.0 # Young's modulus
|
|
variable nu equal 0.3 # Poisson ratio
|
|
variable rho equal 1 # initial mass density
|
|
variable q1 equal 0.56 # standard artificial viscosity linear coefficient
|
|
variable q2 equal 0.0 # standard artificial viscosity quadratic coefficient
|
|
variable hg equal 10.0 # hourglass control coefficient
|
|
variable cp equal 1.0 # heat capacity of material -- not used here
|
|
|
|
variable JC_A equal 0.3241 # Johnson Cook arameters
|
|
variable JC_B equal 0.1138
|
|
variable JC_N equal 0.42
|
|
variable JC_C equal 0 #0.002
|
|
variable JC_M equal 1.34
|
|
variable JC_epsdot0 equal 1.0e-3 # 1/s = 1/(1000 ms)
|
|
variable Troom equal 298.15
|
|
variable Tmelt equal 500.15
|
|
|
|
variable eosC0 equal 0.0 # Polynomial EOS parameters
|
|
variable eosC1 equal 74.2
|
|
variable eosC2 equal 60.5
|
|
variable eosC3 equal 36.5
|
|
variable eosC4 equal 1.96
|
|
variable eosC5 equal 0.0
|
|
variable eosC6 equal 0.0
|
|
|
|
####################################################################################################
|
|
# INITIALIZE LAMMPS
|
|
####################################################################################################
|
|
dimension 2
|
|
units si
|
|
boundary sm sm p # simulation box boundaries
|
|
atom_style smd
|
|
atom_modify map array
|
|
comm_modify vel yes
|
|
neigh_modify every 10 delay 0 check yes # re-build neighbor list every 10 steps
|
|
newton off
|
|
|
|
####################################################################################################
|
|
# CREATE INITIAL GEOMETRY
|
|
####################################################################################################
|
|
variable l0 equal 1.0 # lattice spacing for creating particles
|
|
lattice sq ${l0}
|
|
region box block -10 10 -10 10 -0.1 0.1 units box
|
|
create_box 1 box
|
|
create_atoms 1 box
|
|
group tlsph type 1
|
|
|
|
|
|
####################################################################################################
|
|
# DISCRETIZATION PARAMETERS
|
|
####################################################################################################
|
|
variable h equal 2.01*${l0} # SPH smoothing kernel radius
|
|
variable vol_one equal ${l0}^2 # volume of one particle -- assuming unit thickness
|
|
variable skin equal ${h} # Verlet list range
|
|
neighbor ${skin} bin
|
|
set group all volume ${vol_one}
|
|
set group all smd/mass/density ${rho}
|
|
set group all diameter ${h} # set SPH kernel radius
|
|
|
|
####################################################################################################
|
|
# DEFINE VELOCITY BOUNDARY CONDITIONS
|
|
####################################################################################################
|
|
variable vel0 equal 0.02 # pull velocity
|
|
region top block EDGE EDGE 9.0 EDGE EDGE EDGE units box
|
|
region bot block EDGE EDGE EDGE -9.1 EDGE EDGE units box
|
|
group top region top
|
|
group bot region bot
|
|
variable vel_up equal ${vel0}*(1.0-exp(-0.01*time))
|
|
variable vel_down equal -v_vel_up
|
|
fix veltop_fix top smd/setvelocity 0 v_vel_up 0
|
|
fix velbot_fix bot smd/setvelocity 0 v_vel_down 0
|
|
|
|
####################################################################################################
|
|
# INTERACTION PHYSICS / MATERIAL MODEL
|
|
# We use polynomial EOS for the pressure and the Johnson Cook strength model
|
|
# An integration point fails (cannot support tension anymore) if the plastic strain exceeds 0.5.
|
|
####################################################################################################
|
|
pair_style smd/tlsph
|
|
pair_coeff 1 1 *COMMON ${rho} ${E} ${nu} ${q1} ${q2} ${hg} ${cp} &
|
|
*EOS_POLYNOMIAL ${eosC0} ${eosC1} ${eosC2} ${eosC3} ${eosC4} ${eosC5} ${eosC6} &
|
|
*JOHNSON_COOK ${JC_A} ${JC_B} ${JC_N} ${JC_C} ${JC_epsdot0} ${Troom} ${Tmelt} ${JC_M} &
|
|
*FAILURE_MAX_PLASTIC_STRAIN 1.2 &
|
|
*END
|
|
|
|
####################################################################################################
|
|
# TIME INTEGRATION
|
|
####################################################################################################
|
|
fix dtfix tlsph smd/adjust_dt 0.1 # dynamically adjust time increment every step
|
|
fix integration_fix tlsph smd/integrate_tlsph
|
|
|
|
####################################################################################################
|
|
# SPECIFY TRAJECTORY OUTPUT
|
|
####################################################################################################
|
|
compute dt_atom all smd/tlsph/dt
|
|
compute p all smd/plastic/strain
|
|
compute epsdot all smd/plastic/strain/rate
|
|
compute S all smd/tlsph/stress # Cauchy stress tensor
|
|
compute D all smd/tlsph/strain/rate
|
|
compute E all smd/tlsph/strain
|
|
compute nn all smd/tlsph/num/neighs # number of neighbors for each particle
|
|
compute shape all smd/tlsph/shape
|
|
compute damage all smd/damage
|
|
dump dump_id all custom 100 dump.LAMMPS id type x y z &
|
|
c_S[1] c_S[2] c_S[3] c_S[4] c_S[5] c_S[6] c_S[7] c_nn c_p &
|
|
c_E[1] c_E[2] c_E[3] c_E[4] c_E[5] c_E[6] &
|
|
c_shape[1] c_shape[2] c_shape[3] c_shape[4] c_shape[5] c_shape[6] c_shape[7] &
|
|
c_D[1] c_D[2] c_D[4] c_damage radius c_epsdot &
|
|
vx vy vz c_dt_atom
|
|
dump_modify dump_id first yes
|
|
|
|
####################################################################################################
|
|
# STATUS OUTPUT
|
|
####################################################################################################
|
|
variable stress equal 0.5*(f_velbot_fix[2]-f_veltop_fix[2])/20
|
|
variable length equal xcm(top,y)-xcm(bot,y)
|
|
variable strain equal (v_length-${length})/${length} # engineering strain
|
|
variable time equal f_dtfix
|
|
fix stress_curve all print 10 "${time} ${strain} ${stress}" file stress_strain.dat screen no
|
|
|
|
thermo 100
|
|
thermo_style custom step dt f_dtfix time v_strain
|
|
|
|
####################################################################################################
|
|
# RUN SIMULATION
|
|
####################################################################################################
|
|
#fix 2 all enforce2d
|
|
run 25000
|