Squashed version of Will's commits.
This commit is contained in:
Sachith Dunatunga
46
examples/granular/README_MDR.md
Normal file
46
examples/granular/README_MDR.md
Normal file
@ -0,0 +1,46 @@
|
||||
# Building the executable:
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
CMAKE may be used to build the executable, below is a short shell script that builds a fresh version of the executable.
|
||||
|
||||
#!/bin/bash
|
||||
|
||||
rm -r build
|
||||
mkdir build; cd build
|
||||
cmake ../cmake
|
||||
PKGS="-D PKG_GRANULAR=on -D PKG_VTK=on -D PKG_SRD=on -D PKG_MPI=on"
|
||||
cmake -C ../cmake/presets/most.cmake -C ../cmake/presets/nolib.cmake $PKGS ../cmake
|
||||
cmake --build . -- -j 10
|
||||
|
||||
The GRANULAR package allows DEM simulation in LAMMPS and the VTK package allows creation of VTK files for viewing in paraview. The SRD package is for fluid coupling and the MPI is for parallelization. There are various other ways to build the lammps executable if you are interested: https://docs.lammps.org/Build.html.
|
||||
|
||||
|
||||
# Running the code:
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
In the sims directory of the code base there is a folder called avicelTableting. Inside there is an input file called in.avicelTableting200. This input file defines a tableting simulation that is composed of die filling, compaction, release, and ejection. The geometry of the simulation is modeled after a real compaction simulator die. The material properties that are set roughly represent Avicel PH102. To run this simulation in serial, navigate inside the avicelTableting folder and run the following command:
|
||||
|
||||
(path to lmp executable) < in.avicelTableting200
|
||||
|
||||
The simulation consists of 200 particles and takes ~10 minutes to run on a 2021 M1 Max. To generate a new packing with particles of different radii or distribution you can use the matlab script within generatePackings called generateCylinderPacking.m. It allows you to define the cylinder size, number of desired particles, minimum radius, and maximum radius. It will then generate a packing called spheres.data that is saved in the avicelTableting folder. If you change the distribution of the particle packing make sure to update the following items in the input script:
|
||||
|
||||
neighbor - should be ~1.5x greater than the max radius Rmax.
|
||||
timestep - dt = 0.35*sqrt{m/k} = 0.35*sqrt{(rho*4/3*pi*Rmin^3)/(kappa*Rmin)}.
|
||||
variable atomRadius - specified value in front of prefactor should match Rmin from generateCylinderPacking.m.
|
||||
|
||||
Once running the program should output three files a dump file and two csv files. The dump file can be used in Ovito for visualization. If you desire vtk files for visualization in Paraview create a folder called post and uncomment the dump command in in.avicelTableting. The generated csv files contain the upper punch displacement force relation and particle stress information.
|
||||
|
||||
Plots of the axial stress and radial stress versus displacement can be generated using the matlab script avicelTabletingPlotStresses.m. The script will work even while the simulation is in the middle of running and is a good way to check on the status of the simulation.
|
||||
|
||||
The other input script in the folder is in.avicelTableting20000 this 20,000 particle simulation is the exact one presented in Zunker et al., 2024 and is most efficiently run in parallel with the following command:
|
||||
|
||||
mpirun --np 4 (path to lmp executable) -in in.avicelTableting20000
|
||||
|
||||
Anticipate this simulation to take longer to run, on a 2021 M1 Max running in parallel on 4 processors the simulation took 28 hours to complete.
|
||||
|
||||
Also within the sims directory is the MPFEM folder, which contains two small simulations involving the compaction of 14 monodisperse particles and 12 tridisperse particles as described in Zunker et al., 2024. The input files for these simulations can then be run using the following commands:
|
||||
|
||||
(path to lmp executable) < in.triaxial14particles
|
||||
(path to lmp executable) < in.triaxial12particles
|
||||
|
||||
The outputs of each of these simulations are a dump file and a .csv containing the measured forces in each principal direction. To visualize the measured stresses you can run the Matlab script macro_stresses.m. This will plot stresses for both the LAMMPS results and data from corresponding FEM simulations. The Abaqus 2022 input files for the monodisperse and tridisperse FEM are included in the folder as MPFEM.inp and MPFEM_diff_radii.inp, respectively. For convenience, the principal stress output data for these FEM simulations is already included in the folder as fem.mat and fem_diff_radii.mat.
|
||||
184
examples/granular/in.avicelTableting200
Normal file
184
examples/granular/in.avicelTableting200
Normal file
@ -0,0 +1,184 @@
|
||||
############################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere 1
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 1.0e-3 bin
|
||||
neigh_modify every 10 delay 2000 check no
|
||||
timestep 1.35e-7
|
||||
#processors 2 2 1
|
||||
|
||||
######################### SIMULATION BOUNDING BOX, INTEGRATION, AND, GRAVITY ###########################
|
||||
|
||||
boundary f f f
|
||||
read_data spheres200.data
|
||||
|
||||
######################################### ADD DIE AND ATOM PARAMETERIZATION ##############################################
|
||||
|
||||
variable atomRadius equal 0.44e-3*1.25
|
||||
variable atomDiameter equal 2*${atomRadius}
|
||||
variable atomDensity equal 1560
|
||||
variable atomMassAvg equal ${atomDensity}*4.0/3.0*PI*${atomRadius}^3.0
|
||||
variable dieRadius equal 4e-3
|
||||
variable dieHeight equal 1e-2
|
||||
|
||||
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################
|
||||
|
||||
pair_style granular
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
# linear_history = k_t, x_gamma,t, mu_s
|
||||
variable YoungsModulus equal 5e9
|
||||
variable YieldStress equal 1.9e8
|
||||
variable PoissonsRatio equal 0.4
|
||||
variable SurfaceEnergy equal 2000
|
||||
variable SurfaceEnergyWall equal 0.0
|
||||
variable CoR equal 0.5
|
||||
variable psi_b equal 0.5
|
||||
variable beta equal -ln(${CoR})/sqrt(ln(${CoR})^2+PI^2)
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable kt_wall equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0 #${beta}*sqrt(${atomMassAvg}*${YoungsModulus}*${atomRadius})
|
||||
variable mu_s equal 0.7
|
||||
variable mu_s_wall equal 0.1
|
||||
variable mu_roll equal 0.6
|
||||
variable k_roll equal 2.25*${mu_roll}*${mu_roll}*${YoungsModulus}*${atomRadius}
|
||||
variable gamma_roll equal 0.0 #${beta}*sqrt(${atomMassAvg}*${YoungsModulus}*${atomRadius})
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} damping none tangential linear_history ${kt} ${xgammat} ${mu_s} rolling sds ${k_roll} ${gamma_roll} ${mu_roll}
|
||||
#pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none
|
||||
#pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_nohistory 0.0 0.0 damping none
|
||||
|
||||
######################################### ADD DIE AND PUNCH WALLS ################################################
|
||||
|
||||
variable disp_upper equal 0.0
|
||||
variable disp_lower equal 0.0
|
||||
variable disp_xPlanePos equal 0.0
|
||||
variable disp_xPlaneNeg equal 0.0
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} damping none tangential linear_history ${kt_wall} ${xgammat} ${mu_s_wall} rolling sds ${k_roll} ${gamma_roll} ${mu_roll}"
|
||||
#variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none"
|
||||
#variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} tangential linear_nohistory 0.0 0.0 damping none"
|
||||
|
||||
variable dieHeight2 equal 2*${dieHeight}
|
||||
|
||||
region lowerPunch plane 0 0 0 0 0 1 side in units box move NULL NULL v_disp_lower units box
|
||||
region upperPunch plane 0 0 ${dieHeight} 0 0 -1 side in move NULL NULL v_disp_upper units box
|
||||
region die cylinder z 0 0 ${dieRadius} 0 ${dieHeight2} side in units box
|
||||
variable dieRadiusPlus equal ${dieRadius}*1.05
|
||||
variable dieRadiusPlusNeg equal -${dieRadiusPlus}
|
||||
region xPlanePos plane ${dieRadiusPlus} 0 0 -1 0 0 side in units box move v_disp_xPlanePos NULL NULL units box
|
||||
region xPlaneNeg plane ${dieRadiusPlusNeg} 0 0 1 0 0 side in units box move v_disp_xPlaneNeg NULL NULL units box
|
||||
|
||||
fix lowerPunch all wall/gran/region ${wall_contact_string} region lowerPunch contacts
|
||||
fix upperPunch all wall/gran/region ${wall_contact_string} region upperPunch contacts
|
||||
fix die all wall/gran/region ${wall_contact_string} region die contacts
|
||||
fix xPlanePos all wall/gran/region ${wall_contact_string} region xPlanePos contacts
|
||||
fix xPlaneNeg all wall/gran/region ${wall_contact_string} region xPlaneNeg contacts
|
||||
|
||||
compute avgPunchForce all reduce sum f_upperPunch[4]
|
||||
variable avgPunchForce equal c_avgPunchForce
|
||||
compute avgLowerPunchForce all reduce sum f_lowerPunch[4]
|
||||
variable avgLowerPunchForce equal c_avgLowerPunchForce
|
||||
compute fractureForce all reduce sum f_xPlanePos[2]
|
||||
variable fractureForce equal c_fractureForce
|
||||
|
||||
fix print1 all print 1 "${disp_upper} ${avgPunchForce} ${avgLowerPunchForce} ${fractureForce}" file upperPunchDispForce.csv screen no
|
||||
fix printSTL all print 1 "${disp_upper} ${disp_lower}" file punch_disp_STL_generation.csv screen no
|
||||
|
||||
##################################### INSERT PARTICLES ####################################################
|
||||
|
||||
fix 1 all nve/sphere
|
||||
fix grav all gravity 9.81 vector 0 0 -1
|
||||
|
||||
######################################## SCREEN OUTPUT ####################################################
|
||||
|
||||
#variable ax atom fx/mass
|
||||
#variable ay atom fy/mass
|
||||
#variable az atom fz/mass
|
||||
#variable aMag atom sqrt(v_ax*v_ax+v_ay*v_ay+v_az*v_az)
|
||||
#variable force atom sqrt(fx^2+fy^2+fz^2)
|
||||
#compute aMagAvg all reduce ave v_aMag
|
||||
#variable forceToAccelRatio atom v_force/(mass*v_aMag)
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke vol v_disp_upper #c_aMagAvg
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
compute ke all ke/atom
|
||||
variable output_rate equal round(1e-3/dt)
|
||||
#variable output_rate equal 10
|
||||
|
||||
#dump 1 all custom 1 uniaxialCompression id type f_plane_yz_neg[*]
|
||||
#variable az_upperPunch atom f_upperPunch[7]
|
||||
#variable afz_upperPunch atom f_upperPunch[4]
|
||||
#fix log all print 1 "${delta} ${fx_neg}" file uniaxialCompression.csv title "del,fx_yz_neg" screen no
|
||||
|
||||
run 0
|
||||
|
||||
compute sigmaxx all property/atom d_sigmaxx
|
||||
compute sigmayy all property/atom d_sigmayy
|
||||
compute sigmazz all property/atom d_sigmazz
|
||||
compute Velas all property/atom d_Velas
|
||||
compute adhesive_length all property/atom d_adhesive_length
|
||||
|
||||
compute sigmaxx_ave all reduce ave c_sigmaxx
|
||||
compute sigmayy_ave all reduce ave c_sigmayy
|
||||
compute sigmazz_ave all reduce ave c_sigmazz
|
||||
compute Velas_sum all reduce sum c_Velas
|
||||
compute adhesive_length_ave all reduce ave c_adhesive_length
|
||||
|
||||
variable sxx_ave equal c_sigmaxx_ave
|
||||
variable syy_ave equal c_sigmayy_ave
|
||||
variable szz_ave equal c_sigmazz_ave
|
||||
variable Vparticles equal c_Velas_sum
|
||||
variable adh_length_ave equal c_adhesive_length_ave
|
||||
|
||||
fix log all print 1 "${sxx_ave} ${syy_ave} ${szz_ave} ${Vparticles} ${adh_length_ave}" file avgStresses.csv screen no
|
||||
dump dumpParticles all custom ${output_rate} avicelTableting.dump id type mass diameter x y z vx vy vz fx fy fz c_ke c_sigmaxx c_sigmayy c_sigmazz #v_aMag v_force v_forceToAccelRatio
|
||||
#dump dumpParticlesVTK all vtk ${output_rate} post/particles_*.vtk id x y z fx fy fz vx vy vz c_ke radius c_sigmaxx c_sigmayy c_sigmazz
|
||||
|
||||
######################################### RUN SIMULATION ##########################################
|
||||
|
||||
variable upper_punch_stroke equal 0.6733*${dieHeight} # 0.6733
|
||||
variable vel_upper equal 0.25
|
||||
|
||||
variable settling_steps equal round(0.02/dt)
|
||||
variable compression_steps equal 2*round(${upper_punch_stroke}/${vel_upper}/dt)
|
||||
variable ejection_steps equal ${compression_steps}
|
||||
variable free_float_steps equal round(0.02/dt)
|
||||
variable total_steps equal ${settling_steps}+${compression_steps}+${ejection_steps}+${free_float_steps}
|
||||
|
||||
print "Total steps = ${total_steps}"
|
||||
|
||||
##### SETTLING #####
|
||||
|
||||
run ${settling_steps}
|
||||
|
||||
##### Compression & Release #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${compression_steps}/2)
|
||||
variable disp_upper equal -${upper_punch_stroke}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable short_release equal round(${compression_steps}*1.0)
|
||||
run ${short_release} # 189170 upto #
|
||||
|
||||
##### EJECTION #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${ejection_steps})
|
||||
variable disp_lower equal ${dieHeight}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable disp_upper equal 0.9*v_disp_lower
|
||||
run ${ejection_steps}
|
||||
|
||||
##### FREE FLOAT #####
|
||||
|
||||
variable disp_lower equal ${dieHeight}
|
||||
variable disp_upper equal ${dieHeight}*0.9
|
||||
variable max_disp equal ${dieRadius}*0.75
|
||||
#variable disp_xPlanePos equal -${max_disp}*elapsed/${free_float_steps}
|
||||
#variable disp_xPlaneNeg equal -v_disp_xPlanePos
|
||||
run ${free_float_steps}
|
||||
|
||||
184
examples/granular/in.avicelTableting20000
Normal file
184
examples/granular/in.avicelTableting20000
Normal file
@ -0,0 +1,184 @@
|
||||
############################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere 1
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 0.22e-3 bin
|
||||
neigh_modify every 10 delay 2000 check no
|
||||
timestep 0.3e-7
|
||||
processors 2 2 1
|
||||
|
||||
######################### SIMULATION BOUNDING BOX, INTEGRATION, AND, GRAVITY ###########################
|
||||
|
||||
boundary f f f
|
||||
read_data spheres.data
|
||||
|
||||
######################################### ADD DIE AND ATOM PARAMETERIZATION ##############################################
|
||||
|
||||
variable atomRadius equal 0.0961e-3*1.25
|
||||
variable atomDiameter equal 2*${atomRadius}
|
||||
variable atomDensity equal 1560
|
||||
variable atomMassAvg equal ${atomDensity}*4.0/3.0*PI*${atomRadius}^3.0
|
||||
variable dieRadius equal 4e-3
|
||||
variable dieHeight equal 1e-2
|
||||
|
||||
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################
|
||||
|
||||
pair_style granular
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
# linear_history = k_t, x_gamma,t, mu_s
|
||||
variable YoungsModulus equal 5e9
|
||||
variable YieldStress equal 1.9e8
|
||||
variable PoissonsRatio equal 0.4
|
||||
variable SurfaceEnergy equal 450
|
||||
variable SurfaceEnergyWall equal 0.0
|
||||
variable CoR equal 0.5
|
||||
variable psi_b equal 0.5
|
||||
variable beta equal -ln(${CoR})/sqrt(ln(${CoR})^2+PI^2)
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable kt_wall equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0 #${beta}*sqrt(${atomMassAvg}*${YoungsModulus}*${atomRadius})
|
||||
variable mu_s equal 0.7
|
||||
variable mu_s_wall equal 0.1
|
||||
variable mu_roll equal 0.6
|
||||
variable k_roll equal 2.25*${mu_roll}*${mu_roll}*${YoungsModulus}*${atomRadius}
|
||||
variable gamma_roll equal 0.0 #${beta}*sqrt(${atomMassAvg}*${YoungsModulus}*${atomRadius})
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} damping none tangential linear_history ${kt} ${xgammat} ${mu_s} rolling sds ${k_roll} ${gamma_roll} ${mu_roll}
|
||||
#pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none
|
||||
#pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_nohistory 0.0 0.0 damping none
|
||||
|
||||
######################################### ADD DIE AND PUNCH WALLS ################################################
|
||||
|
||||
variable disp_upper equal 0.0
|
||||
variable disp_lower equal 0.0
|
||||
variable disp_xPlanePos equal 0.0
|
||||
variable disp_xPlaneNeg equal 0.0
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} damping none tangential linear_history ${kt_wall} ${xgammat} ${mu_s_wall} rolling sds ${k_roll} ${gamma_roll} ${mu_roll}"
|
||||
#variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none"
|
||||
#variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} tangential linear_nohistory 0.0 0.0 damping none"
|
||||
|
||||
variable dieHeight2 equal 2*${dieHeight}
|
||||
|
||||
region lowerPunch plane 0 0 0 0 0 1 side in units box move NULL NULL v_disp_lower units box
|
||||
region upperPunch plane 0 0 ${dieHeight} 0 0 -1 side in move NULL NULL v_disp_upper units box
|
||||
region die cylinder z 0 0 ${dieRadius} 0 ${dieHeight2} side in units box
|
||||
variable dieRadiusPlus equal ${dieRadius}*1.05
|
||||
variable dieRadiusPlusNeg equal -${dieRadiusPlus}
|
||||
region xPlanePos plane ${dieRadiusPlus} 0 0 -1 0 0 side in units box move v_disp_xPlanePos NULL NULL units box
|
||||
region xPlaneNeg plane ${dieRadiusPlusNeg} 0 0 1 0 0 side in units box move v_disp_xPlaneNeg NULL NULL units box
|
||||
|
||||
fix lowerPunch all wall/gran/region ${wall_contact_string} region lowerPunch contacts
|
||||
fix upperPunch all wall/gran/region ${wall_contact_string} region upperPunch contacts
|
||||
fix die all wall/gran/region ${wall_contact_string} region die contacts
|
||||
fix xPlanePos all wall/gran/region ${wall_contact_string} region xPlanePos contacts
|
||||
fix xPlaneNeg all wall/gran/region ${wall_contact_string} region xPlaneNeg contacts
|
||||
|
||||
compute avgPunchForce all reduce sum f_upperPunch[4]
|
||||
variable avgPunchForce equal c_avgPunchForce
|
||||
compute avgLowerPunchForce all reduce sum f_lowerPunch[4]
|
||||
variable avgLowerPunchForce equal c_avgLowerPunchForce
|
||||
compute fractureForce all reduce sum f_xPlanePos[2]
|
||||
variable fractureForce equal c_fractureForce
|
||||
|
||||
fix print1 all print 1 "${disp_upper} ${avgPunchForce} ${avgLowerPunchForce} ${fractureForce}" file upperPunchDispForce.csv screen no
|
||||
fix printSTL all print 1 "${disp_upper} ${disp_lower}" file punch_disp_STL_generation.csv screen no
|
||||
|
||||
##################################### INSERT PARTICLES ####################################################
|
||||
|
||||
fix 1 all nve/sphere
|
||||
fix grav all gravity 9.81 vector 0 0 -1
|
||||
|
||||
######################################## SCREEN OUTPUT ####################################################
|
||||
|
||||
#variable ax atom fx/mass
|
||||
#variable ay atom fy/mass
|
||||
#variable az atom fz/mass
|
||||
#variable aMag atom sqrt(v_ax*v_ax+v_ay*v_ay+v_az*v_az)
|
||||
#variable force atom sqrt(fx^2+fy^2+fz^2)
|
||||
#compute aMagAvg all reduce ave v_aMag
|
||||
#variable forceToAccelRatio atom v_force/(mass*v_aMag)
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke vol v_disp_upper #c_aMagAvg
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
compute ke all ke/atom
|
||||
variable output_rate equal round(1e-4/dt)
|
||||
#variable output_rate equal 10
|
||||
|
||||
#dump 1 all custom 1 uniaxialCompression id type f_plane_yz_neg[*]
|
||||
#variable az_upperPunch atom f_upperPunch[7]
|
||||
#variable afz_upperPunch atom f_upperPunch[4]
|
||||
#fix log all print 1 "${delta} ${fx_neg}" file uniaxialCompression.csv title "del,fx_yz_neg" screen no
|
||||
|
||||
run 0
|
||||
|
||||
compute sigmaxx all property/atom d_sigmaxx
|
||||
compute sigmayy all property/atom d_sigmayy
|
||||
compute sigmazz all property/atom d_sigmazz
|
||||
compute Velas all property/atom d_Velas
|
||||
compute adhesive_length all property/atom d_adhesive_length
|
||||
|
||||
compute sigmaxx_ave all reduce ave c_sigmaxx
|
||||
compute sigmayy_ave all reduce ave c_sigmayy
|
||||
compute sigmazz_ave all reduce ave c_sigmazz
|
||||
compute Velas_sum all reduce sum c_Velas
|
||||
compute adhesive_length_ave all reduce ave c_adhesive_length
|
||||
|
||||
variable sxx_ave equal c_sigmaxx_ave
|
||||
variable syy_ave equal c_sigmayy_ave
|
||||
variable szz_ave equal c_sigmazz_ave
|
||||
variable Vparticles equal c_Velas_sum
|
||||
variable adh_length_ave equal c_adhesive_length_ave
|
||||
|
||||
fix log all print 1 "${sxx_ave} ${syy_ave} ${szz_ave} ${Vparticles} ${adh_length_ave}" file avgStresses.csv screen no
|
||||
dump dumpParticles all custom ${output_rate} avicelTableting.dump id type mass diameter x y z vx vy vz fx fy fz c_ke c_sigmaxx c_sigmayy c_sigmazz #v_aMag v_force v_forceToAccelRatio
|
||||
#dump dumpParticlesVTK all vtk ${output_rate} post/particles_*.vtk id x y z fx fy fz vx vy vz c_ke radius c_sigmaxx c_sigmayy c_sigmazz
|
||||
|
||||
######################################### RUN SIMULATION ##########################################
|
||||
|
||||
variable upper_punch_stroke equal 0.6905*${dieHeight} # 0.6733
|
||||
variable vel_upper equal 0.25
|
||||
|
||||
variable settling_steps equal round(0.02/dt)
|
||||
variable compression_steps equal 2*round(${upper_punch_stroke}/${vel_upper}/dt)
|
||||
variable ejection_steps equal ${compression_steps}
|
||||
variable free_float_steps equal round(0.02/dt)
|
||||
variable total_steps equal ${settling_steps}+${compression_steps}+${ejection_steps}+${free_float_steps}
|
||||
|
||||
print "Total steps = ${total_steps}"
|
||||
|
||||
##### SETTLING #####
|
||||
|
||||
run ${settling_steps}
|
||||
|
||||
##### Compression & Release #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${compression_steps}/2)
|
||||
variable disp_upper equal -${upper_punch_stroke}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable short_release equal round(${compression_steps}*1.0)
|
||||
run ${short_release} # 189170 upto #
|
||||
|
||||
##### EJECTION #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${ejection_steps})
|
||||
variable disp_lower equal ${dieHeight}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable disp_upper equal 0.9*v_disp_lower
|
||||
run ${ejection_steps}
|
||||
|
||||
##### FREE FLOAT #####
|
||||
|
||||
variable disp_lower equal ${dieHeight}
|
||||
variable disp_upper equal ${dieHeight}*0.9
|
||||
variable max_disp equal ${dieRadius}*0.75
|
||||
#variable disp_xPlanePos equal -${max_disp}*elapsed/${free_float_steps}
|
||||
#variable disp_xPlaneNeg equal -v_disp_xPlanePos
|
||||
run ${free_float_steps}
|
||||
|
||||
100
examples/granular/in.mpfem.triaxial12particles
Normal file
100
examples/granular/in.mpfem.triaxial12particles
Normal file
@ -0,0 +1,100 @@
|
||||
############################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 2 bin
|
||||
neigh_modify delay 0
|
||||
timestep 1e-6
|
||||
|
||||
######################### SIMULATION BOUNDING BOX, INTEGRATION, AND, GRAVITY ###########################
|
||||
|
||||
boundary f f f
|
||||
read_data spheres.data
|
||||
fix integr all nve/sphere
|
||||
|
||||
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################
|
||||
|
||||
variable atomRadius equal 0.5
|
||||
|
||||
pair_style granular
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
# linear_history = k_t, x_gamma,t, mu_s
|
||||
variable YoungsModulus equal 1e9
|
||||
variable PoissonsRatio equal 0.3
|
||||
variable YieldStress equal 50e6
|
||||
variable SurfaceEnergy equal 0.0
|
||||
variable psi_b equal 0.5
|
||||
variable CoR equal 0.5
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0
|
||||
variable mu_s equal 0.5
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none
|
||||
|
||||
######################################### ADD IN PLANES ################################################
|
||||
|
||||
variable boxWidth equal 3
|
||||
variable halfBoxWidth equal ${boxWidth}/2
|
||||
|
||||
variable plane_disp equal 0.0
|
||||
variable plane_disp_neg equal 0.0
|
||||
|
||||
region plane_yz_pos plane ${halfBoxWidth} 0 0 -1 0 0 side in move v_plane_disp_neg NULL NULL units box
|
||||
region plane_yz_neg plane -${halfBoxWidth} 0 0 1 0 0 side in move v_plane_disp NULL NULL units box
|
||||
region plane_xz_pos plane 0 ${halfBoxWidth} 0 0 -1 0 side in move NULL v_plane_disp_neg NULL units box
|
||||
region plane_xz_neg plane 0 -${halfBoxWidth} 0 0 1 0 side in move NULL v_plane_disp NULL units box
|
||||
region plane_xy_pos plane 0 0 ${halfBoxWidth} 0 0 -1 side in move NULL NULL v_plane_disp_neg units box
|
||||
region plane_xy_neg plane 0 0 -${halfBoxWidth} 0 0 1 side in move NULL NULL v_plane_disp units box
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none"
|
||||
|
||||
fix plane_yz_pos all wall/gran/region ${wall_contact_string} region plane_yz_pos contacts
|
||||
fix plane_yz_neg all wall/gran/region ${wall_contact_string} region plane_yz_neg contacts
|
||||
fix plane_xz_pos all wall/gran/region ${wall_contact_string} region plane_xz_pos contacts
|
||||
fix plane_xz_neg all wall/gran/region ${wall_contact_string} region plane_xz_neg contacts
|
||||
fix plane_xy_pos all wall/gran/region ${wall_contact_string} region plane_xy_pos contacts
|
||||
fix plane_xy_neg all wall/gran/region ${wall_contact_string} region plane_xy_neg contacts
|
||||
|
||||
compute plane_xy_neg_force all reduce sum f_plane_xy_neg[4]
|
||||
variable plane_xy_neg_force equal c_plane_xy_neg_force
|
||||
|
||||
compute plane_xz_neg_force all reduce sum f_plane_xz_neg[3]
|
||||
variable plane_xz_neg_force equal c_plane_xz_neg_force
|
||||
|
||||
compute plane_yz_neg_force all reduce sum f_plane_yz_neg[2]
|
||||
variable plane_yz_neg_force equal c_plane_yz_neg_force
|
||||
|
||||
fix print1 all print 1 "${plane_disp} ${plane_xy_neg_force} ${plane_xz_neg_force} ${plane_yz_neg_force}" file lammps_macro_force_disp_diff_radii.csv screen no
|
||||
|
||||
######################################## SCREEN OUTPUT ####################################################
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke c_1 vol
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### DEFINE WALL MOVEMENT #################################################
|
||||
|
||||
variable disp_max equal 0.499
|
||||
variable ddisp equal 0.00001
|
||||
variable compression_steps equal round(${disp_max}/${ddisp})
|
||||
variable output_rate equal round(${compression_steps}/100)
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
#dump dmp all vtk ${output_rate} post/triaxial12particles_*.vtk id type mass x y z vx vy vz fx fy fz radius
|
||||
dump dumpParticles all custom ${output_rate} triaxial12particles.dump id type mass x y z vx vy vz fx fy fz radius
|
||||
|
||||
#################################### COMPRESS THE PARTICLES ##################################################
|
||||
|
||||
run 0
|
||||
|
||||
variable plane_disp equal ${ddisp}*elapsed
|
||||
variable plane_disp_neg equal -${ddisp}*elapsed
|
||||
|
||||
run ${compression_steps}
|
||||
|
||||
#################################### UNCOMPRESS THE PARTICLES ################################################
|
||||
164
examples/granular/in.mpfem.triaxial14particles
Normal file
164
examples/granular/in.mpfem.triaxial14particles
Normal file
@ -0,0 +1,164 @@
|
||||
############################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere 1
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 2 bin
|
||||
neigh_modify delay 0
|
||||
timestep 1e-6
|
||||
#processors 1 2 1
|
||||
|
||||
######################### SIMULATION BOUNDING BOX, INTEGRATION, AND, GRAVITY ###########################
|
||||
|
||||
boundary f f f
|
||||
region box block -100 100 -100 100 -100 100 units box
|
||||
create_box 2 box
|
||||
fix integr all nve/sphere
|
||||
|
||||
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################
|
||||
|
||||
variable atomRadius equal 0.5
|
||||
|
||||
pair_style granular
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
# linear_history = k_t, x_gamma,t, mu_s
|
||||
variable YoungsModulus equal 1e9
|
||||
variable PoissonsRatio equal 0.3
|
||||
variable YieldStress equal 50e6
|
||||
variable SurfaceEnergy equal 0.0
|
||||
variable psi_b equal 0.5
|
||||
variable CoR equal 0.5
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0
|
||||
variable mu_s equal 0.5
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none
|
||||
|
||||
######################################### ADD IN PLANES ################################################
|
||||
|
||||
variable boxWidth equal 3
|
||||
variable halfBoxWidth equal ${boxWidth}/2
|
||||
|
||||
variable plane_disp equal 0.0
|
||||
variable plane_disp_neg equal 0.0
|
||||
|
||||
region plane_yz_pos plane ${halfBoxWidth} 0 0 -1 0 0 side in move v_plane_disp_neg NULL NULL units box
|
||||
region plane_yz_neg plane -${halfBoxWidth} 0 0 1 0 0 side in move v_plane_disp NULL NULL units box
|
||||
region plane_xz_pos plane 0 ${halfBoxWidth} 0 0 -1 0 side in move NULL v_plane_disp_neg NULL units box
|
||||
region plane_xz_neg plane 0 -${halfBoxWidth} 0 0 1 0 side in move NULL v_plane_disp NULL units box
|
||||
region plane_xy_pos plane 0 0 ${halfBoxWidth} 0 0 -1 side in move NULL NULL v_plane_disp_neg units box
|
||||
region plane_xy_neg plane 0 0 -${halfBoxWidth} 0 0 1 side in move NULL NULL v_plane_disp units box
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none"
|
||||
|
||||
fix plane_yz_pos all wall/gran/region ${wall_contact_string} region plane_yz_pos contacts
|
||||
fix plane_yz_neg all wall/gran/region ${wall_contact_string} region plane_yz_neg contacts
|
||||
fix plane_xz_pos all wall/gran/region ${wall_contact_string} region plane_xz_pos contacts
|
||||
fix plane_xz_neg all wall/gran/region ${wall_contact_string} region plane_xz_neg contacts
|
||||
fix plane_xy_pos all wall/gran/region ${wall_contact_string} region plane_xy_pos contacts
|
||||
fix plane_xy_neg all wall/gran/region ${wall_contact_string} region plane_xy_neg contacts
|
||||
|
||||
compute plane_xy_neg_force all reduce sum f_plane_xy_neg[4]
|
||||
variable plane_xy_neg_force equal c_plane_xy_neg_force
|
||||
|
||||
compute plane_xz_neg_force all reduce sum f_plane_xz_neg[3]
|
||||
variable plane_xz_neg_force equal c_plane_xz_neg_force
|
||||
|
||||
compute plane_yz_neg_force all reduce sum f_plane_yz_neg[2]
|
||||
variable plane_yz_neg_force equal c_plane_yz_neg_force
|
||||
|
||||
fix print1 all print 1 "${plane_disp} ${plane_xy_neg_force} ${plane_xz_neg_force} ${plane_yz_neg_force}" file lammps_macro_force_disp.csv screen no
|
||||
|
||||
##################################### INSERT PARTICLES ####################################################
|
||||
|
||||
# PARTICLE PROPERTIES
|
||||
variable d equal 1
|
||||
variable den equal 1000
|
||||
|
||||
create_atoms 1 single -0.1682 0.1113 0.1492 units box
|
||||
group atomOne id 1
|
||||
set atom 1 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.6234 -0.5134 0.4198 units box
|
||||
group atomTwo id 2
|
||||
set atom 2 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.8958 0.8304 0.1353 units box
|
||||
group atomThree id 3
|
||||
set atom 3 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.1058 -0.9321 -0.2803 units box
|
||||
group atomFour id 4
|
||||
set atom 4 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.7816 0.5089 -0.5358 units box
|
||||
group atomFive id 5
|
||||
set atom 5 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.09169 -0.1675 -0.9413 units box
|
||||
group atomSix id 6
|
||||
set atom 6 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.8324 -0.6242 -0.9531 units box
|
||||
group atomSeven id 7
|
||||
set atom 7 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.8994 -0.6587 0.9025 units box
|
||||
group atomEight id 8
|
||||
set atom 8 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.08605 0.8565 0.8629 units box
|
||||
group atomNine id 9
|
||||
set atom 9 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.7989 0.8445 -0.9942 units box
|
||||
group atomTen id 10
|
||||
set atom 10 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.9317 -0.8971 -0.8523 units box
|
||||
group atomEleven id 11
|
||||
set atom 11 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single -0.979 0.3867 0.7977 units box
|
||||
group atomTwelve id 12
|
||||
set atom 12 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.9926 0.2609 0.9693 units box
|
||||
group atomThirteen id 13
|
||||
set atom 13 diameter ${d} density ${den}
|
||||
|
||||
create_atoms 1 single 0.9814 0.01161 -0.4335 units box
|
||||
group atomFourteen id 14
|
||||
set atom 14 diameter ${d} density ${den}
|
||||
|
||||
######################################## SCREEN OUTPUT ####################################################
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke c_1 vol
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### DEFINE WALL MOVEMENT #################################################
|
||||
|
||||
variable disp_max equal 0.45
|
||||
variable ddisp equal 0.00001
|
||||
variable compression_steps equal round(${disp_max}/${ddisp})
|
||||
variable output_rate equal round(${compression_steps}/100)
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
#dump dmp all vtk ${output_rate} post/triaxial14particles_*.vtk id type mass x y z vx vy vz fx fy fz radius
|
||||
dump dumpParticles all custom ${output_rate} triaxial14particles.dump id type mass x y z vx vy vz fx fy fz radius
|
||||
|
||||
#################################### COMPRESS THE PARTICLES ##################################################
|
||||
|
||||
run 0
|
||||
|
||||
variable plane_disp equal ${ddisp}*elapsed
|
||||
variable plane_disp_neg equal -${ddisp}*elapsed
|
||||
|
||||
run ${compression_steps}
|
||||
|
||||
#################################### UNCOMPRESS THE PARTICLES ################################################
|
||||
23
examples/granular/spheres.data
Normal file
23
examples/granular/spheres.data
Normal file
@ -0,0 +1,23 @@
|
||||
#LAMMPS data file created by matlab.
|
||||
12 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-10.0000000000 10.0000000000 xlo xhi
|
||||
-10.0000000000 10.0000000000 ylo yhi
|
||||
-10.0000000000 10.0000000000 zlo zhi
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 0.8000000000 1000.0000000000 0.0717535226 -0.2092222842 0.3662146798
|
||||
2 1 1.2000000000 1000.0000000000 -0.8233763986 -0.7426114800 -0.8263932264
|
||||
3 1 0.8000000000 1000.0000000000 -1.0685863278 -0.4494609702 0.2196698078
|
||||
4 1 0.8000000000 1000.0000000000 0.5829432471 -1.0098803839 -0.7607543861
|
||||
5 1 0.8000000000 1000.0000000000 -0.8658471132 0.6951192569 0.0107556658
|
||||
6 1 1.2000000000 1000.0000000000 0.3966456126 0.7215053869 -0.7540113087
|
||||
7 1 1.2000000000 1000.0000000000 0.7316242921 0.8996483982 0.6751483031
|
||||
8 1 1.0000000000 1000.0000000000 0.6267527768 -0.8419367233 0.6964197101
|
||||
9 1 0.8000000000 1000.0000000000 -0.0409043189 -0.1452314035 -1.0102948313
|
||||
10 1 0.8000000000 1000.0000000000 -0.9495107709 0.6760151650 -0.9220534482
|
||||
11 1 1.0000000000 1000.0000000000 -0.7488486472 0.2188003421 0.7892021020
|
||||
12 1 1.2000000000 1000.0000000000 0.8968590780 -0.2350366437 -0.2006719701
|
||||
211
examples/granular/spheres200.data
Normal file
211
examples/granular/spheres200.data
Normal file
@ -0,0 +1,211 @@
|
||||
#LAMMPS data file created by matlab.
|
||||
200 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-0.005000 0.005000 xlo xhi
|
||||
-0.005000 0.005000 ylo yhi
|
||||
-0.001000 0.020000 zlo zhi
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 0.001206 1560.000000 -0.000938 0.000556 0.000883
|
||||
2 1 0.000953 1560.000000 -0.002626 -0.000145 0.002778
|
||||
3 1 0.001035 1560.000000 -0.000434 0.000172 0.008458
|
||||
4 1 0.001225 1560.000000 -0.003126 -0.000604 0.004986
|
||||
5 1 0.001119 1560.000000 0.000772 0.002972 0.002568
|
||||
6 1 0.001243 1560.000000 -0.000363 0.001184 0.004927
|
||||
7 1 0.001173 1560.000000 0.000218 0.000243 0.005475
|
||||
8 1 0.000937 1560.000000 0.000033 0.000029 0.003141
|
||||
9 1 0.001055 1560.000000 -0.001660 0.001975 0.008611
|
||||
10 1 0.000938 1560.000000 -0.001818 0.002352 0.002534
|
||||
11 1 0.000990 1560.000000 0.001592 0.000435 0.004416
|
||||
12 1 0.000927 1560.000000 -0.001659 -0.000004 0.005901
|
||||
13 1 0.001272 1560.000000 0.002972 0.000553 0.007291
|
||||
14 1 0.001226 1560.000000 0.002090 0.000983 0.001406
|
||||
15 1 0.000957 1560.000000 0.002241 -0.001608 0.001304
|
||||
16 1 0.001020 1560.000000 -0.001944 0.001290 0.002030
|
||||
17 1 0.001289 1560.000000 -0.002256 -0.001173 0.003474
|
||||
18 1 0.000998 1560.000000 0.000771 0.002127 0.000906
|
||||
19 1 0.000927 1560.000000 0.000186 0.000567 0.001207
|
||||
20 1 0.001095 1560.000000 -0.000937 -0.003179 0.008173
|
||||
21 1 0.001006 1560.000000 -0.001736 0.000751 0.004618
|
||||
22 1 0.001037 1560.000000 0.000784 0.001844 0.002380
|
||||
23 1 0.001297 1560.000000 0.000234 -0.001597 0.008560
|
||||
24 1 0.001017 1560.000000 0.002454 -0.000505 0.001171
|
||||
25 1 0.001110 1560.000000 -0.000803 -0.000415 0.003714
|
||||
26 1 0.001192 1560.000000 0.002283 0.000648 0.003048
|
||||
27 1 0.000992 1560.000000 -0.000065 -0.000545 0.007062
|
||||
28 1 0.001116 1560.000000 0.002174 -0.001463 0.005830
|
||||
29 1 0.001258 1560.000000 0.001602 0.001853 0.007246
|
||||
30 1 0.001055 1560.000000 -0.001535 -0.002770 0.007196
|
||||
31 1 0.000958 1560.000000 -0.000438 -0.000260 0.004709
|
||||
32 1 0.001188 1560.000000 0.000339 -0.000355 0.009171
|
||||
33 1 0.001166 1560.000000 0.002513 -0.001215 0.004434
|
||||
34 1 0.000907 1560.000000 0.001905 -0.000373 0.004921
|
||||
35 1 0.001245 1560.000000 -0.000091 -0.002620 0.004150
|
||||
36 1 0.001302 1560.000000 0.003292 0.000184 0.005377
|
||||
37 1 0.001305 1560.000000 0.002099 0.001261 0.008939
|
||||
38 1 0.000988 1560.000000 0.003274 0.000136 0.003667
|
||||
39 1 0.000892 1560.000000 0.001798 -0.002104 0.008610
|
||||
40 1 0.001247 1560.000000 -0.003058 -0.000575 0.000948
|
||||
41 1 0.000900 1560.000000 -0.000258 -0.000469 0.001478
|
||||
42 1 0.000945 1560.000000 -0.001434 -0.001711 0.004610
|
||||
43 1 0.000977 1560.000000 -0.001410 0.002808 0.004963
|
||||
44 1 0.000930 1560.000000 -0.002110 -0.001362 0.006749
|
||||
45 1 0.000931 1560.000000 0.001256 -0.000876 0.000844
|
||||
46 1 0.000901 1560.000000 0.000899 -0.001189 0.005316
|
||||
47 1 0.000940 1560.000000 -0.002189 -0.000047 0.007240
|
||||
48 1 0.001217 1560.000000 -0.000108 -0.001333 0.002257
|
||||
49 1 0.001088 1560.000000 0.001364 -0.000594 0.002789
|
||||
50 1 0.001143 1560.000000 -0.000311 -0.001425 0.006092
|
||||
51 1 0.001054 1560.000000 0.002262 0.002312 0.004315
|
||||
52 1 0.001016 1560.000000 -0.000724 0.000741 0.003295
|
||||
53 1 0.001051 1560.000000 0.000527 -0.001987 0.003307
|
||||
54 1 0.000905 1560.000000 0.000827 0.001457 0.005868
|
||||
55 1 0.001195 1560.000000 -0.001176 -0.000645 0.000798
|
||||
56 1 0.001253 1560.000000 0.002583 -0.001847 0.003310
|
||||
57 1 0.000982 1560.000000 0.001551 -0.002803 0.005076
|
||||
58 1 0.000945 1560.000000 -0.000481 0.000354 0.007220
|
||||
59 1 0.001040 1560.000000 -0.002736 0.001076 0.008769
|
||||
60 1 0.000917 1560.000000 0.000826 -0.001887 0.006449
|
||||
61 1 0.000914 1560.000000 -0.001171 -0.001592 0.007266
|
||||
62 1 0.000959 1560.000000 0.000834 -0.002671 0.007105
|
||||
63 1 0.000990 1560.000000 -0.000251 -0.001327 0.004339
|
||||
64 1 0.001220 1560.000000 0.001384 0.002896 0.005874
|
||||
65 1 0.000949 1560.000000 -0.001340 -0.000608 0.007496
|
||||
66 1 0.001306 1560.000000 0.002187 0.002068 0.002629
|
||||
67 1 0.001206 1560.000000 0.000148 0.001506 0.008517
|
||||
68 1 0.001123 1560.000000 0.001288 -0.000303 0.006613
|
||||
69 1 0.001151 1560.000000 -0.000876 0.001549 0.001740
|
||||
70 1 0.001315 1560.000000 -0.001902 -0.002590 0.001344
|
||||
71 1 0.000927 1560.000000 0.002285 -0.000866 0.006900
|
||||
72 1 0.001279 1560.000000 -0.000165 0.002689 0.007449
|
||||
73 1 0.000910 1560.000000 0.001009 0.001054 0.005049
|
||||
74 1 0.001148 1560.000000 -0.002229 -0.001285 0.008736
|
||||
75 1 0.001067 1560.000000 -0.000261 -0.002945 0.002157
|
||||
76 1 0.000993 1560.000000 -0.001641 0.002272 0.007601
|
||||
77 1 0.001228 1560.000000 0.001939 -0.000214 0.008903
|
||||
78 1 0.001076 1560.000000 0.000767 0.001172 0.003556
|
||||
79 1 0.001105 1560.000000 -0.000561 0.002493 0.004214
|
||||
80 1 0.001195 1560.000000 0.002694 -0.000817 0.007949
|
||||
81 1 0.001239 1560.000000 -0.000968 -0.003145 0.006096
|
||||
82 1 0.001083 1560.000000 -0.000808 0.001813 0.006396
|
||||
83 1 0.000923 1560.000000 0.000632 -0.001437 0.001310
|
||||
84 1 0.000981 1560.000000 -0.001842 0.002774 0.006508
|
||||
85 1 0.000998 1560.000000 -0.002775 0.001616 0.001453
|
||||
86 1 0.000979 1560.000000 -0.002520 0.001715 0.007741
|
||||
87 1 0.001002 1560.000000 -0.001465 -0.001931 0.006048
|
||||
88 1 0.000958 1560.000000 0.003264 0.000707 0.001189
|
||||
89 1 0.001052 1560.000000 -0.001314 -0.000701 0.002721
|
||||
90 1 0.001096 1560.000000 0.001154 0.002129 0.004403
|
||||
91 1 0.001104 1560.000000 0.002118 0.001977 0.000794
|
||||
92 1 0.001263 1560.000000 -0.001499 -0.002764 0.003441
|
||||
93 1 0.001086 1560.000000 -0.001096 0.002514 0.001154
|
||||
94 1 0.000895 1560.000000 0.001130 0.000029 0.001045
|
||||
95 1 0.000964 1560.000000 0.000905 -0.003200 0.000542
|
||||
96 1 0.000898 1560.000000 -0.000868 0.003148 0.008306
|
||||
97 1 0.000907 1560.000000 -0.001406 0.001144 0.007862
|
||||
98 1 0.001176 1560.000000 0.001246 -0.001074 0.004327
|
||||
99 1 0.001148 1560.000000 0.001512 -0.002739 0.003346
|
||||
100 1 0.000922 1560.000000 0.001470 -0.000036 0.007695
|
||||
101 1 0.001031 1560.000000 -0.002751 0.000928 0.004124
|
||||
102 1 0.001030 1560.000000 -0.000177 -0.002370 0.005374
|
||||
103 1 0.000915 1560.000000 0.000824 0.000521 0.007070
|
||||
104 1 0.001085 1560.000000 -0.002281 -0.000023 0.009123
|
||||
105 1 0.001004 1560.000000 -0.000167 0.002610 0.008905
|
||||
106 1 0.001060 1560.000000 -0.000389 -0.002220 0.007688
|
||||
107 1 0.000920 1560.000000 -0.000483 0.003231 0.006505
|
||||
108 1 0.001122 1560.000000 0.001781 -0.001547 0.002237
|
||||
109 1 0.001172 1560.000000 -0.002650 0.000830 0.005429
|
||||
110 1 0.001137 1560.000000 -0.000030 -0.003246 0.001024
|
||||
111 1 0.001315 1560.000000 0.001470 -0.001735 0.007580
|
||||
112 1 0.001245 1560.000000 0.000481 -0.003067 0.006025
|
||||
113 1 0.000904 1560.000000 0.000632 -0.000184 0.002010
|
||||
114 1 0.000883 1560.000000 -0.001828 0.002191 0.003819
|
||||
115 1 0.000974 1560.000000 0.002167 0.001616 0.006226
|
||||
116 1 0.001150 1560.000000 0.000871 -0.002731 0.002136
|
||||
117 1 0.001312 1560.000000 -0.000326 -0.001971 0.001000
|
||||
118 1 0.000914 1560.000000 0.001020 0.000810 0.002086
|
||||
119 1 0.001136 1560.000000 -0.000101 -0.003277 0.007246
|
||||
120 1 0.000991 1560.000000 -0.001944 0.000576 0.003215
|
||||
121 1 0.001216 1560.000000 -0.000913 -0.001165 0.008857
|
||||
122 1 0.001045 1560.000000 -0.003110 0.001062 0.002973
|
||||
123 1 0.000918 1560.000000 0.000348 0.000365 0.004046
|
||||
124 1 0.001279 1560.000000 -0.000884 0.003087 0.002268
|
||||
125 1 0.001065 1560.000000 -0.002238 0.001309 0.006452
|
||||
126 1 0.001012 1560.000000 -0.002059 -0.001354 0.001935
|
||||
127 1 0.001142 1560.000000 -0.003011 0.000567 0.001739
|
||||
128 1 0.000921 1560.000000 0.001764 0.002804 0.008177
|
||||
129 1 0.001151 1560.000000 -0.003105 -0.000384 0.006602
|
||||
130 1 0.000967 1560.000000 0.000932 0.000588 0.008823
|
||||
131 1 0.000908 1560.000000 -0.001873 -0.001947 0.007825
|
||||
132 1 0.000923 1560.000000 -0.002993 0.000883 0.007425
|
||||
133 1 0.001171 1560.000000 0.003310 -0.000405 0.006558
|
||||
134 1 0.000977 1560.000000 -0.000098 -0.000180 0.000492
|
||||
135 1 0.000938 1560.000000 -0.000706 -0.000129 0.006085
|
||||
136 1 0.001008 1560.000000 -0.000256 0.002333 0.000550
|
||||
137 1 0.001073 1560.000000 0.000534 -0.000055 0.008080
|
||||
138 1 0.000890 1560.000000 0.000351 0.001695 0.007195
|
||||
139 1 0.000973 1560.000000 0.002593 0.001907 0.005394
|
||||
140 1 0.001176 1560.000000 -0.001862 -0.000534 0.004494
|
||||
141 1 0.001306 1560.000000 -0.000951 0.001053 0.009299
|
||||
142 1 0.001103 1560.000000 -0.001937 -0.002711 0.008485
|
||||
143 1 0.001262 1560.000000 -0.002947 -0.001470 0.007682
|
||||
144 1 0.000914 1560.000000 0.002047 0.000811 0.005504
|
||||
145 1 0.000954 1560.000000 0.001935 -0.002349 0.006632
|
||||
146 1 0.001003 1560.000000 0.000766 -0.002635 0.008483
|
||||
147 1 0.001137 1560.000000 0.000102 0.003195 0.004922
|
||||
148 1 0.001006 1560.000000 -0.001982 0.001014 0.000685
|
||||
149 1 0.001255 1560.000000 -0.000718 0.001939 0.003056
|
||||
150 1 0.001057 1560.000000 -0.001189 -0.001717 0.003045
|
||||
151 1 0.001228 1560.000000 0.001581 0.002926 0.003510
|
||||
152 1 0.001052 1560.000000 -0.002172 0.001949 0.004831
|
||||
153 1 0.000979 1560.000000 -0.001817 0.000291 0.002048
|
||||
154 1 0.001286 1560.000000 -0.002647 -0.001839 0.004620
|
||||
155 1 0.001085 1560.000000 -0.000081 0.000850 0.002139
|
||||
156 1 0.000990 1560.000000 -0.000081 0.002105 0.005587
|
||||
157 1 0.001043 1560.000000 0.001636 -0.000112 0.001860
|
||||
158 1 0.001309 1560.000000 0.003216 -0.000851 0.002791
|
||||
159 1 0.000913 1560.000000 0.000608 0.003148 0.006565
|
||||
160 1 0.000919 1560.000000 0.000536 -0.003106 0.003249
|
||||
161 1 0.000943 1560.000000 0.003145 -0.000528 0.008915
|
||||
162 1 0.000993 1560.000000 -0.002811 -0.000099 0.008110
|
||||
163 1 0.001125 1560.000000 0.001415 -0.002271 0.000643
|
||||
164 1 0.000919 1560.000000 -0.001406 0.000223 0.006781
|
||||
165 1 0.001040 1560.000000 0.000690 0.003193 0.008329
|
||||
166 1 0.001055 1560.000000 0.001075 0.002584 0.009093
|
||||
167 1 0.001176 1560.000000 0.000851 0.003176 0.000591
|
||||
168 1 0.001003 1560.000000 -0.001462 0.001511 0.005544
|
||||
169 1 0.001126 1560.000000 -0.000077 0.003324 0.001347
|
||||
170 1 0.001068 1560.000000 0.003110 0.000810 0.008495
|
||||
171 1 0.001011 1560.000000 -0.001661 0.000117 0.008201
|
||||
172 1 0.001066 1560.000000 -0.000359 -0.003279 0.009094
|
||||
173 1 0.001303 1560.000000 0.003066 0.001188 0.004082
|
||||
174 1 0.000983 1560.000000 0.000354 0.002261 0.003558
|
||||
175 1 0.001137 1560.000000 0.002860 -0.001571 0.009180
|
||||
176 1 0.001070 1560.000000 0.001246 -0.001279 0.009104
|
||||
177 1 0.000886 1560.000000 0.002271 -0.000316 0.003675
|
||||
178 1 0.000983 1560.000000 -0.001987 -0.002490 0.005377
|
||||
179 1 0.000939 1560.000000 0.000601 -0.000861 0.003477
|
||||
180 1 0.001177 1560.000000 0.001522 0.002902 0.001690
|
||||
181 1 0.001036 1560.000000 -0.001200 -0.002874 0.004750
|
||||
182 1 0.000898 1560.000000 -0.001705 -0.001140 0.005503
|
||||
183 1 0.001315 1560.000000 0.002732 0.001766 0.007885
|
||||
184 1 0.001318 1560.000000 -0.002909 -0.001610 0.005936
|
||||
185 1 0.001218 1560.000000 0.003213 0.000884 0.002316
|
||||
186 1 0.001234 1560.000000 -0.002394 -0.002298 0.002575
|
||||
187 1 0.001160 1560.000000 -0.003313 -0.000065 0.003625
|
||||
188 1 0.001022 1560.000000 -0.003096 -0.001048 0.002151
|
||||
189 1 0.000966 1560.000000 0.001891 -0.002093 0.004404
|
||||
190 1 0.001048 1560.000000 -0.002367 0.002338 0.000697
|
||||
191 1 0.000995 1560.000000 -0.001204 -0.001912 0.002030
|
||||
192 1 0.001136 1560.000000 -0.001152 -0.002402 0.009223
|
||||
193 1 0.001083 1560.000000 -0.002588 -0.001768 0.000753
|
||||
194 1 0.000946 1560.000000 -0.001338 -0.000741 0.006527
|
||||
195 1 0.000943 1560.000000 -0.000073 0.003254 0.003663
|
||||
196 1 0.001059 1560.000000 0.000087 0.000958 0.006388
|
||||
197 1 0.001131 1560.000000 0.001030 0.001019 0.000752
|
||||
198 1 0.001257 1560.000000 -0.001365 0.002946 0.009266
|
||||
199 1 0.000891 1560.000000 -0.000445 -0.000273 0.002382
|
||||
200 1 0.001055 1560.000000 0.001781 0.000748 0.006583
|
||||
20011
examples/granular/spheres20000.data
Normal file
20011
examples/granular/spheres20000.data
Normal file
File diff suppressed because it is too large
Load Diff
21
examples/granular/spheresSTLgeneration.data
Normal file
21
examples/granular/spheresSTLgeneration.data
Normal file
@ -0,0 +1,21 @@
|
||||
#LAMMPS data file created by matlab.
|
||||
10 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-0.005000 0.005000 xlo xhi
|
||||
-0.005000 0.005000 ylo yhi
|
||||
-0.001000 0.020000 zlo zhi
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 0.000249 1560.000000 0.000469 -0.000120 0.008380
|
||||
2 1 0.000260 1560.000000 -0.001811 0.000425 0.004090
|
||||
3 1 0.000225 1560.000000 -0.001371 0.000426 0.002952
|
||||
4 1 0.000252 1560.000000 0.000928 -0.000089 0.004118
|
||||
5 1 0.000210 1560.000000 0.000208 0.000384 0.000842
|
||||
6 1 0.000218 1560.000000 -0.002456 -0.000587 0.003834
|
||||
7 1 0.000201 1560.000000 -0.000278 0.003459 0.003445
|
||||
8 1 0.000231 1560.000000 -0.002123 0.003188 0.004799
|
||||
9 1 0.000268 1560.000000 -0.000338 -0.001557 0.007382
|
||||
10 1 0.000244 1560.000000 0.002537 0.002151 0.005999
|
||||
@ -496,6 +496,7 @@ void FixWallGran::post_force(int /*vflag*/)
|
||||
model->dx[2] = dz;
|
||||
model->radi = radius[i];
|
||||
model->radj = rwall;
|
||||
|
||||
if (model->beyond_contact) model->touch = history_one[i][0];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
|
||||
@ -50,14 +50,15 @@ class FixWallGran : public Fix {
|
||||
int maxsize_restart() override;
|
||||
void reset_dt() override;
|
||||
|
||||
// for granular model choices
|
||||
class Granular_NS::GranularModel *model; // MOVED HERE FROM PROTECTED FOR MDR MODEL
|
||||
void clear_stored_contacts(); // MOVED HERE FROM PROTECTED FOR MDR MODEL
|
||||
|
||||
protected:
|
||||
int wallstyle, wiggle, wshear, axis;
|
||||
int nlevels_respa;
|
||||
bigint time_origin;
|
||||
|
||||
// for granular model choices
|
||||
class Granular_NS::GranularModel *model;
|
||||
|
||||
double lo, hi, cylradius;
|
||||
double amplitude, period, omega, vshear;
|
||||
double dt;
|
||||
@ -86,7 +87,7 @@ class FixWallGran : public Fix {
|
||||
|
||||
int store;
|
||||
|
||||
void clear_stored_contacts();
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -29,6 +29,7 @@
|
||||
#include "region.h"
|
||||
#include "update.h"
|
||||
#include "variable.h"
|
||||
#include <iostream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
@ -129,6 +130,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
if (update->setupflag) history_update = 0;
|
||||
model->history_update = history_update;
|
||||
|
||||
//std::cout << (uint64_t)this << ", " << (uint64_t)model << std::endl;
|
||||
|
||||
// if just reneighbored:
|
||||
// update rigid body masses for owned atoms if using FixRigid
|
||||
// body[i] = which body atom I is in, -1 if none
|
||||
@ -228,6 +231,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
model->radi = radius[i];
|
||||
model->radj = region->contact[ic].radius;
|
||||
model->r = region->contact[ic].r;
|
||||
model->i = i; // Added for MDR
|
||||
model->j = ic; // Added for MDR
|
||||
|
||||
if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
|
||||
@ -44,23 +44,29 @@ class FixWallGranRegion : public FixWallGran {
|
||||
int size_restart(int) override;
|
||||
int maxsize_restart() override;
|
||||
|
||||
class Region *region; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
void update_contacts(int, int); // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
int *ncontact; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
int **walls; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
int *c2r; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
double ***history_many; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
int tmax; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
|
||||
|
||||
private:
|
||||
class Region *region;
|
||||
|
||||
int nregion;
|
||||
|
||||
// shear history for multiple contacts per particle
|
||||
|
||||
int tmax; // max # of region walls one particle can touch
|
||||
int *ncontact; // # of shear contacts per particle
|
||||
int **walls; // which wall each contact is with
|
||||
double ***history_many; // history per particle per contact
|
||||
int *c2r; // contact to region mapping
|
||||
|
||||
|
||||
|
||||
// c2r[i] = index of Ith contact in
|
||||
// region-contact[] list of contacts
|
||||
int motion_resetflag; // used by restart to indicate that region
|
||||
// vel info is to be reset
|
||||
|
||||
void update_contacts(int, int);
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -16,8 +16,14 @@
|
||||
#include "error.h"
|
||||
#include "granular_model.h"
|
||||
#include "math_const.h"
|
||||
#include "atom.h"
|
||||
#include "csv_writer.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <sstream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
@ -381,3 +387,767 @@ void GranSubModNormalJKR::set_fncrit()
|
||||
{
|
||||
Fncrit = fabs(Fne + 2.0 * F_pulloff);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
MDR contact model
|
||||
|
||||
Contributing authors:
|
||||
William Zunker (MIT), Sachith Dunatunga (MIT),
|
||||
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
|
||||
GranSubModNormal(gm, lmp)
|
||||
{
|
||||
num_coeffs = 6; // Young's Modulus, Poisson's ratio, yield stress, effective surface energy, psi_b, coefficent of restitution
|
||||
contact_radius_flag = 1;
|
||||
size_history = 26;
|
||||
|
||||
nondefault_history_transfer = 1;
|
||||
transfer_history_factor = new double[size_history];
|
||||
//transfer_history_factor[0] = +1;
|
||||
for (int i = 0; i < size_history; i++) {
|
||||
transfer_history_factor[i] = +1;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void GranSubModNormalMDR::coeffs_to_local()
|
||||
{
|
||||
E = coeffs[0]; // Young's modulus
|
||||
nu = coeffs[1]; // Poisson's ratio
|
||||
Y = coeffs[2]; // yield stress
|
||||
gamma = coeffs[3]; // effective surface energy
|
||||
psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total}
|
||||
CoR = coeffs[5]; // coefficent of restitution
|
||||
|
||||
if (E <= 0.0) error->all(FLERR, "Illegal MDR normal model, Young's modulus must be greater than 0");
|
||||
if (nu < 0.0 || nu > 0.5) error->all(FLERR, "Illegal MDR normal model, Poisson's ratio must be between 0 and 0.5");
|
||||
if (Y < 0.0) error->all(FLERR, "Illegal MDR normal model, yield stress must be greater than or equal to 0");
|
||||
if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0");
|
||||
if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0");
|
||||
if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double GranSubModNormalMDR::calculate_forces()
|
||||
|
||||
{
|
||||
// To understand the structure of the overall code it is important to consider
|
||||
// the following:
|
||||
//
|
||||
// The MDR contact model was developed by imagining individual particles being
|
||||
// squished between a number of rigid flats (references below). To allow
|
||||
// for many interacting particles, we extend the idea of isolated particles surrounded
|
||||
// by rigid flats. In particular, we imagine placing rigid flats at the overlap
|
||||
// midpoints between particles. The force is calculated seperately on both sides
|
||||
// of the contact assuming interaction with a rigid flat. The two forces are then
|
||||
// averaged on either side of the contact to determine the final force. If the
|
||||
// contact is between a particle and wall then only one force evaluation is required.
|
||||
//
|
||||
// Zunker and Kamrin, 2024, Part I: https://doi.org/10.1016/j.jmps.2023.105492
|
||||
// Zunker and Kamrin, 2024, Part II: https://doi.org/10.1016/j.jmps.2023.105493
|
||||
|
||||
const int itag_true = atom->tag[gm->i]; // true i particle tag
|
||||
const int jtag_true = atom->tag[gm->j]; // true j particle tag
|
||||
const int i_true = gm->i; // true i particle index
|
||||
const int j_true = gm->j; // true j particle index
|
||||
const double radi_true = gm->radi; // true i particle initial radius
|
||||
const double radj_true = gm->radj; // true j particle initial radius
|
||||
|
||||
F = 0.0; // average force
|
||||
double F0 = 0.0; // force on contact side 0
|
||||
double F1 = 0.0; // force on contact side 1
|
||||
double R0 = 0.0;
|
||||
double R1 = 0.0;
|
||||
int i0 = 0;
|
||||
int i1 = 0;
|
||||
double delta = gm->delta; // apparent overlap
|
||||
|
||||
//if (gm->contact_type == PAIR) delta = gm->delta/2.0; // half displacement to imagine interaction with rigid flat
|
||||
//std::cout << "Normal force is called for: " << i_true << ", " << j_true << std::endl;
|
||||
//std::cout << "Contact model has been entered " << gm->contact_type << ", " << PAIR << ", " << WALL << ", " << WALLREGION << ", " << gm->itype << ", " << gm->jtype << ", " << gm->delta << std::endl;
|
||||
|
||||
// initialize indexing in history array of different constact history variables
|
||||
const int delta_offset_0 = 0; // apparent overlap
|
||||
const int delta_offset_1 = 1;
|
||||
const int deltao_offset_0 = 2; // displacement
|
||||
const int deltao_offset_1 = 3;
|
||||
const int delta_MDR_offset_0 = 4; // MDR apparent overlap
|
||||
const int delta_MDR_offset_1 = 5;
|
||||
const int delta_BULK_offset_0 = 6; // bulk displacement
|
||||
const int delta_BULK_offset_1 = 7;
|
||||
const int deltamax_MDR_offset_0 = 8; // maximum MDR apparent overlap
|
||||
const int deltamax_MDR_offset_1 = 9;
|
||||
const int Yflag_offset_0 = 10; // yield flag
|
||||
const int Yflag_offset_1 = 11;
|
||||
const int deltaY_offset_0 = 12; // yield displacement
|
||||
const int deltaY_offset_1 = 13;
|
||||
const int cA_offset_0 = 14; // contact area intercept
|
||||
const int cA_offset_1 = 15;
|
||||
const int aAdh_offset_0 = 16; // adhesive contact radius
|
||||
const int aAdh_offset_1 = 17;
|
||||
const int Ac_offset_0 = 18; // contact area
|
||||
const int Ac_offset_1 = 19;
|
||||
const int eps_bar_offset_0 = 20; // volume-averaged infinitesimal strain tensor
|
||||
const int eps_bar_offset_1 = 21;
|
||||
const int penalty_offset_ = 22; // contact penalty
|
||||
const int deltamax_offset_ = 23;
|
||||
const int deltap_offset_0 = 24;
|
||||
const int deltap_offset_1 = 25;
|
||||
|
||||
|
||||
// initialize particle history variables
|
||||
int tmp1, tmp2;
|
||||
int index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
|
||||
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
|
||||
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
|
||||
int index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
|
||||
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
|
||||
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
|
||||
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
|
||||
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
|
||||
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
|
||||
int index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
|
||||
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
|
||||
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
|
||||
int index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
|
||||
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation
|
||||
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
|
||||
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
|
||||
int index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
|
||||
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
|
||||
double * Rinitial = atom->dvector[index_Ro];
|
||||
double * Vgeo = atom->dvector[index_Vgeo];
|
||||
double * Velas = atom->dvector[index_Velas];
|
||||
double * Vcaps = atom->dvector[index_Vcaps];
|
||||
double * eps_bar = atom->dvector[index_eps_bar];
|
||||
double * dRnumerator = atom->dvector[index_dRnumerator];
|
||||
double * dRdenominator = atom->dvector[index_dRdenominator];
|
||||
double * Acon0 = atom->dvector[index_Acon0];
|
||||
double * Acon1 = atom->dvector[index_Acon1];
|
||||
double * Atot = atom->dvector[index_Atot];
|
||||
double * Atot_sum = atom->dvector[index_Atot_sum];
|
||||
double * ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
double * psi = atom->dvector[index_psi];
|
||||
double * sigmaxx = atom->dvector[index_sigmaxx];
|
||||
double * sigmayy = atom->dvector[index_sigmayy];
|
||||
double * sigmazz = atom->dvector[index_sigmazz];
|
||||
double * contacts = atom->dvector[index_contacts];
|
||||
double * adhesive_length = atom->dvector[index_adhesive_length];
|
||||
|
||||
|
||||
double * history = & gm->history[history_index]; // load in all history variables
|
||||
|
||||
// Rigid flat placement scheme
|
||||
double * deltamax_offset = & history[deltamax_offset_];
|
||||
double deltamax = *deltamax_offset;
|
||||
double * deltap_offset0 = & history[deltap_offset_0];
|
||||
double * deltap_offset1 = & history[deltap_offset_1];
|
||||
double deltap0 = *deltap_offset0;
|
||||
double deltap1 = *deltap_offset1;
|
||||
|
||||
if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta;
|
||||
deltamax = *deltamax_offset;
|
||||
|
||||
for (int contactSide = 0; contactSide < 2; contactSide++) {
|
||||
|
||||
double * delta_offset;
|
||||
double * deltao_offset;
|
||||
double * delta_MDR_offset;
|
||||
double * delta_BULK_offset;
|
||||
double * deltamax_MDR_offset;
|
||||
double * Yflag_offset;
|
||||
double * deltaY_offset;
|
||||
double * cA_offset;
|
||||
double * aAdh_offset;
|
||||
double * Ac_offset;
|
||||
double * eps_bar_offset;
|
||||
double * penalty_offset;
|
||||
double * deltap_offset;
|
||||
|
||||
double overlap_limit = 0.75;
|
||||
|
||||
if (contactSide == 0) {
|
||||
if (gm->contact_type == PAIR) {
|
||||
if (itag_true > jtag_true) {
|
||||
gm->i = i_true;
|
||||
gm->j = j_true;
|
||||
gm->radi = radi_true;
|
||||
gm->radj = radj_true;
|
||||
} else {
|
||||
gm->i = j_true;
|
||||
gm->j = i_true;
|
||||
gm->radi = radj_true;
|
||||
gm->radj = radi_true;
|
||||
}
|
||||
R0 = gm->radi;
|
||||
R1 = gm->radj;
|
||||
i0 = gm->i;
|
||||
i1 = gm->j;
|
||||
|
||||
double delta_geo;
|
||||
double delta_geo_alt;
|
||||
double delta_geoOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1));
|
||||
double delta_geoOpt2 = deltamax*(deltamax - 2.0*R0)/(2.0*(deltamax - R0 - R1));
|
||||
(gm->radi < gm->radj) ? delta_geo = MAX(delta_geoOpt1,delta_geoOpt2) : delta_geo = MIN(delta_geoOpt1,delta_geoOpt2);
|
||||
(gm->radi > gm->radj) ? delta_geo_alt = MAX(delta_geoOpt1,delta_geoOpt2) : delta_geo_alt = MIN(delta_geoOpt1,delta_geoOpt2);
|
||||
|
||||
if (delta_geo/gm->radi > overlap_limit) {
|
||||
delta_geo = gm->radi*overlap_limit;
|
||||
} else if (delta_geo_alt/gm->radj > overlap_limit) {
|
||||
delta_geo = deltamax - gm->radj*overlap_limit;
|
||||
}
|
||||
|
||||
double deltap = deltap0 + deltap1;
|
||||
delta = delta_geo + (deltap0 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
|
||||
|
||||
//std::cout << "CS 0: " << gm->radi << ", " << gm->radj << ", gm->delta " << gm->delta << ", delta " << delta << ", deltamax " << deltamax << ", delta_geo " << delta_geo << ", delta_geo_alt " << delta_geo_alt << ", delta_geo/Ri " << delta_geo/gm->radi << ", delta_geo_alt/Rj " << delta_geo_alt/gm->radj << ", delta_geoOpt1 " << delta_geoOpt1 << ", delta_geoOpt2 " << delta_geoOpt2 << ", deltamax-delta " << (deltamax-gm->delta) << std::endl;
|
||||
}
|
||||
delta_offset = & history[delta_offset_0];
|
||||
deltao_offset = & history[deltao_offset_0];
|
||||
delta_MDR_offset = & history[delta_MDR_offset_0];
|
||||
delta_BULK_offset = & history[delta_BULK_offset_0];
|
||||
deltamax_MDR_offset = & history[deltamax_MDR_offset_0];
|
||||
Yflag_offset = & history[Yflag_offset_0];
|
||||
deltaY_offset = & history[deltaY_offset_0];
|
||||
cA_offset = & history[cA_offset_0];
|
||||
aAdh_offset = & history[aAdh_offset_0];
|
||||
Ac_offset = & history[Ac_offset_0];
|
||||
eps_bar_offset = & history[eps_bar_offset_0];
|
||||
deltap_offset = & history[deltap_offset_0];
|
||||
} else {
|
||||
if (gm->contact_type != PAIR) break; // contact with particle-wall requires only one evaluation
|
||||
if (itag_true < jtag_true) {
|
||||
gm->i = i_true;
|
||||
gm->j = j_true;
|
||||
gm->radi = radi_true;
|
||||
gm->radj = radj_true;
|
||||
} else {
|
||||
gm->i = j_true;
|
||||
gm->j = i_true;
|
||||
gm->radi = radj_true;
|
||||
gm->radj = radi_true;
|
||||
}
|
||||
|
||||
double delta_geo;
|
||||
double delta_geo_alt;
|
||||
double delta_geoOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1));
|
||||
double delta_geoOpt2 = deltamax*(deltamax - 2.0*R0)/(2.0*(deltamax - R0 - R1));
|
||||
(gm->radi < gm->radj) ? delta_geo = MAX(delta_geoOpt1,delta_geoOpt2) : delta_geo = MIN(delta_geoOpt1,delta_geoOpt2);
|
||||
(gm->radi > gm->radj) ? delta_geo_alt = MAX(delta_geoOpt1,delta_geoOpt2) : delta_geo_alt = MIN(delta_geoOpt1,delta_geoOpt2);
|
||||
|
||||
if (delta_geo/gm->radi > overlap_limit) {
|
||||
delta_geo = gm->radi*overlap_limit;
|
||||
} else if (delta_geo_alt/gm->radj > overlap_limit) {
|
||||
delta_geo = deltamax - gm->radj*overlap_limit;
|
||||
}
|
||||
|
||||
double deltap = deltap0 + deltap1;
|
||||
delta = delta_geo + (deltap1 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
|
||||
|
||||
//std::cout << "CS 1: " << gm->radi << ", " << gm->radj << ", gm->delta " << gm->delta << ", delta " << delta << ", deltamax " << deltamax << ", delta_geo " << delta_geo << ", delta_geo_alt " << delta_geo_alt << ", delta_geo/Ri " << delta_geo/gm->radi << ", delta_geo_alt/Rj " << delta_geo_alt/gm->radj << ", delta_geoOpt1 " << delta_geoOpt1 << ", delta_geoOpt2 " << delta_geoOpt2 << ", deltamax-delta " << (deltamax-gm->delta) << std::endl;
|
||||
|
||||
delta_offset = & history[delta_offset_1];
|
||||
deltao_offset = & history[deltao_offset_1];
|
||||
delta_MDR_offset = & history[delta_MDR_offset_1];
|
||||
delta_BULK_offset = & history[delta_BULK_offset_1];
|
||||
deltamax_MDR_offset = & history[deltamax_MDR_offset_1];
|
||||
Yflag_offset = & history[Yflag_offset_1];
|
||||
deltaY_offset = & history[deltaY_offset_1];
|
||||
cA_offset = & history[cA_offset_1];
|
||||
aAdh_offset = & history[aAdh_offset_1];
|
||||
Ac_offset = & history[Ac_offset_1];
|
||||
eps_bar_offset = & history[eps_bar_offset_1];
|
||||
deltap_offset = & history[deltap_offset_1];
|
||||
}
|
||||
|
||||
// temporary i and j indices
|
||||
const int i = gm->i;
|
||||
const int j = gm->j;
|
||||
|
||||
//std::cout << lmp->update->ntimestep << std::endl;
|
||||
|
||||
// material and geometric property definitions
|
||||
// E, nu, Y gamma , psi_b, and CoR are already defined.
|
||||
const double G = E/(2.0*(1.0+nu)); // shear modulus
|
||||
const double kappa = E/(3.0*(1.0-2.0*nu)); // bulk modulus
|
||||
const double Eeff = E/(1.0-pow(nu,2.0)); // composite plane strain modulus
|
||||
|
||||
const double Ro = Rinitial[i]; // initial radius
|
||||
const double R = gm->radi; // apparent radius
|
||||
|
||||
// kinematics
|
||||
const double ddelta = delta - *delta_offset;
|
||||
*delta_offset = delta;
|
||||
|
||||
const double deltao = delta - (R - Ro);
|
||||
const double ddeltao = deltao - *deltao_offset;
|
||||
*deltao_offset = deltao;
|
||||
|
||||
double ddelta_MDR;
|
||||
double ddelta_BULK;
|
||||
if ( psi[i] < psi_b ) { // if true, bulk response has triggered, split displacement increment between the MDR and BULK components
|
||||
ddelta_MDR = std::min(ddelta-ddelta_bar[i], delta-*delta_MDR_offset);
|
||||
ddelta_BULK = ddelta_bar[i];
|
||||
} else { // if false, no bulk response, full displacement increment goes to the MDR component
|
||||
ddelta_BULK = 0.0;
|
||||
ddelta_MDR = ddelta;
|
||||
}
|
||||
const double delta_MDR = *delta_MDR_offset + ddelta_MDR; // MDR displacement
|
||||
*delta_MDR_offset = delta_MDR; // Update old MDR displacement
|
||||
const double delta_BULK = std::max(0.0,*delta_BULK_offset+ddelta_BULK); // bulk displacement
|
||||
*delta_BULK_offset = delta_BULK; // update old bulk displacement
|
||||
|
||||
if (delta_MDR > *deltamax_MDR_offset) *deltamax_MDR_offset = delta_MDR;
|
||||
const double deltamax_MDR = *deltamax_MDR_offset;
|
||||
|
||||
const double pY = Y*(1.75*exp(-4.4*deltamax_MDR/R) + 1.0); // Set value of average pressure along yield surface
|
||||
if ( *Yflag_offset == 0.0 && delta_MDR >= deltamax_MDR ) {
|
||||
const double phertz = 4*Eeff*sqrt(delta_MDR)/(3*M_PI*sqrt(R));
|
||||
if ( phertz > pY ) {
|
||||
*Yflag_offset = 1.0;
|
||||
*deltaY_offset = delta_MDR;
|
||||
*cA_offset = M_PI*(pow(*deltaY_offset,2.0) - *deltaY_offset*R);
|
||||
}
|
||||
}
|
||||
|
||||
//if (i_true == 167 && j_true == 204) {
|
||||
//std::cout << "i " << i << " | j " << j << " | delta_BULK: " << delta_BULK << " | delta_MDR " << delta_MDR << " | ddelta_BULK " << ddelta_BULK << " | ddelta_MDR " << ddelta_MDR << std::endl;
|
||||
//}
|
||||
|
||||
//std::cout << "Yield Flag: " << *Yflag_offset << ", " << R << std::endl;
|
||||
|
||||
// MDR force calculation
|
||||
double F_MDR;
|
||||
double A; // height of elliptical indenter
|
||||
double B; // width of elliptical indenter
|
||||
double deltae1D; // transformed elastic displacement
|
||||
double deltaR; // displacement correction
|
||||
double amax; // maximum experienced contact radius
|
||||
const double cA = *cA_offset; // contact area intercept
|
||||
|
||||
if ( *Yflag_offset == 0.0 ) { // elastic contact
|
||||
A = 4.0*R;
|
||||
B = 2.0*R;
|
||||
deltae1D = delta_MDR;
|
||||
(deltae1D > 0) ? amax = sqrt(deltae1D*R) : amax = 0.0;
|
||||
} else { // plastic contact
|
||||
amax = sqrt((2.0*deltamax_MDR*R - pow(deltamax_MDR,2.0)) + cA/M_PI);
|
||||
A = 4.0*pY/Eeff*amax;
|
||||
B = 2.0*amax;
|
||||
const double deltae1Dmax = A/2.0; // maximum transformed elastic displacement
|
||||
const double Fmax = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1Dmax/A) - (1.0 - 2.0*deltae1Dmax/A)*sqrt(4.0*deltae1Dmax/A - 4.0*pow(deltae1Dmax,2.0)/pow(A,2.0))); // force caused by full submersion of elliptical indenter to depth of A/2
|
||||
const double zR = R - (deltamax_MDR - deltae1Dmax); // depth of particle center
|
||||
deltaR = (Fmax*(2*pow(amax,2.0)*(-1 + nu) - (-1 + 2*nu)*zR*(-zR + sqrt(pow(amax,2.0) + pow(zR,2.0)))))/((M_PI*pow(amax,2.0))*2*G*sqrt(pow(amax,2.0) + pow(zR,2.0)));
|
||||
deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR)/(1 + deltaR/deltae1Dmax); // transformed elastic displacement
|
||||
|
||||
// added for rigid flat placement
|
||||
*deltap_offset = deltamax_MDR - (deltae1Dmax + deltaR);
|
||||
//std::cout << *deltap_offset << std::endl;
|
||||
}
|
||||
|
||||
//std::cout << psi_b << ", " << psi[i] << ", " << A << ", " << B << ", " << pY << ", " << amax << " || " << deltao << ", " << delta << ", " << ddelta << ", " << *delta_offset << ", " << ddelta_bar[i] << " || " << delta_MDR << ", " << ddelta_MDR << ", " << *delta_MDR_offset << ", " << deltamax_MDR << " || " << delta_BULK << ", " << ddelta_BULK << ", " << *delta_BULK_offset << " || " << R << std::endl;
|
||||
//std::cout << i << ", " << j << ", " << A << ", " << B << " || " << deltao << ", " << delta << ", " << ddelta << ", " << R << ", " << M_PI*pow(amax,2.0) << std::endl;
|
||||
|
||||
double a_na;
|
||||
double a_fac = 0.99;
|
||||
(deltae1D >= 0.0) ? a_na = B*sqrt(A - deltae1D)*sqrt(deltae1D)/A : a_na = 0.0;
|
||||
double aAdh = *aAdh_offset;
|
||||
if (aAdh > a_fac*amax) aAdh = a_fac*amax;
|
||||
|
||||
//if (i_true == 4 && j_true == 52){
|
||||
//std::cout << "CS: " << contactSide << ", aAdh: " << aAdh << ", deltae1D: " << deltae1D << ", A: " << A << ", B:" << B << ", amax: " << amax << ", deltae1D: " << deltae1D << ", R: " << R << std::endl;
|
||||
//}
|
||||
|
||||
if ( gamma > 0.0 ) { // adhesive contact
|
||||
double g_aAdh;
|
||||
|
||||
if (delta_MDR == deltamax_MDR || a_na >= aAdh ) { // case 1: no tensile springs, purely compressive contact
|
||||
(deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0)));
|
||||
if ( std::isnan(F_MDR) ) {
|
||||
std::cout << "F_MDR is NaN, case 1: no tensile springs" << std::endl;
|
||||
//std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl;
|
||||
std::cout << "i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", contact type: " << gm->contact_type << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << std::endl;
|
||||
std::exit(1);
|
||||
}
|
||||
*aAdh_offset = a_fac*a_na;
|
||||
} else {
|
||||
const double lmax = sqrt(2.0*M_PI*aAdh*gamma/Eeff);
|
||||
g_aAdh = A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh,2.0));
|
||||
const double acrit = (-((pow(B,2)*gamma*M_PI)/(pow(A,2)*Eeff)) + (pow(2,0.3333333333333333)*pow(B,4)*pow(gamma,2)*pow(M_PI,1.6666666666666667))/
|
||||
(pow(A,2)*pow(Eeff,2)*pow((27*pow(A,4)*pow(B,4)*gamma)/Eeff - (2*pow(B,6)*pow(gamma,3)*pow(M_PI,2))/pow(Eeff,3) + (3*sqrt(3)*sqrt(27*pow(A,8)*pow(B,8)*pow(Eeff,2)*pow(gamma,2) -
|
||||
4*pow(A,4)*pow(B,10)*pow(gamma,4)*pow(M_PI,2)))/pow(Eeff,2),0.3333333333333333)) + (pow(M_PI/2.,0.3333333333333333)*pow((27*pow(A,4)*pow(B,4)*gamma)/Eeff -
|
||||
(2*pow(B,6)*pow(gamma,3)*pow(M_PI,2))/pow(Eeff,3) + (3*sqrt(3)*sqrt(27*pow(A,8)*pow(B,8)*pow(Eeff,2)*pow(gamma,2) - 4*pow(A,4)*pow(B,10)*pow(gamma,4)*pow(M_PI,2)))/
|
||||
pow(Eeff,2),0.3333333333333333))/pow(A,2))/6;
|
||||
|
||||
if ( (deltae1D + lmax - g_aAdh) >= 0.0) { // case 2: tensile springs, but not exceeding critical length --> deltae + lmax - g(aAdhes) >= 0
|
||||
//if (contactSide == 0) {
|
||||
//std::cout << "Case 2 tensile springs not exceeding critical length, R " << R << "deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh" << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl;
|
||||
//}
|
||||
const double deltaeAdh = g_aAdh;
|
||||
const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0)));
|
||||
const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh;
|
||||
F_MDR = F_na + F_Adhes;
|
||||
if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 2: tensile springs, but not exceeding critical length" << std::endl;
|
||||
|
||||
} else { // case 3: tensile springs exceed critical length --> deltae + lmax - g(aAdhes) = 0
|
||||
//if (contactSide == 0) {
|
||||
//std::cout << "Case 3 tensile springs exceed critical length, R " << R << " deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh " << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl;
|
||||
//}
|
||||
|
||||
if ( aAdh < acrit ) {
|
||||
aAdh = 0.0;
|
||||
F_MDR = 0.0;
|
||||
} else {
|
||||
// newton-raphson to find aAdh
|
||||
const double maxIterations = 100;
|
||||
const double error = 1e-10;
|
||||
const double error2 = 1e-16;
|
||||
double aAdh_tmp = aAdh;
|
||||
double fa;
|
||||
double fa2;
|
||||
double dfda;
|
||||
for (int lv1 = 0; lv1 < maxIterations; ++lv1) {
|
||||
fa = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) );
|
||||
if (abs(fa) < error) {
|
||||
break;
|
||||
}
|
||||
dfda = -((aAdh_tmp*A)/(B*sqrt(-pow(aAdh_tmp,2.0) + pow(B,2.0)/4.0))) + (gamma*sqrt(M_PI/2.0))/(Eeff*sqrt((aAdh_tmp*gamma)/Eeff));
|
||||
aAdh_tmp = aAdh_tmp - fa/dfda;
|
||||
fa2 = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) );
|
||||
if (abs(fa-fa2) < error2) {
|
||||
break;
|
||||
}
|
||||
if (lv1 == maxIterations-1){
|
||||
aAdh_tmp = 0.0;
|
||||
}
|
||||
}
|
||||
aAdh = aAdh_tmp;
|
||||
|
||||
g_aAdh = A/2.0 - A/B*sqrt(pow(B,2.0)/4.0 - pow(aAdh,2.0));
|
||||
const double deltaeAdh = g_aAdh;
|
||||
const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0)));
|
||||
const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh;
|
||||
F_MDR = F_na + F_Adhes;
|
||||
if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 3: tensile springs exceed critical length" << std::endl;
|
||||
}
|
||||
*aAdh_offset = aAdh;
|
||||
}
|
||||
}
|
||||
} else { // non-adhesive contact
|
||||
*aAdh_offset = a_na;
|
||||
(deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0)));
|
||||
if ( std::isnan(F_MDR) ) {
|
||||
std::cout << "F_MDR is NaN, non-adhesive case" << std::endl;
|
||||
//std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl;
|
||||
std::cout << "i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << std::endl;
|
||||
std::exit(1);
|
||||
}
|
||||
}
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", aAdh_offset: " << *aAdh_offset << ", aAdh: " << aAdh << ", a_na: " << a_na << std::endl;
|
||||
|
||||
contacts[i] += 1;
|
||||
adhesive_length[i] += aAdh;
|
||||
|
||||
// contact penalty scheme
|
||||
penalty_offset = & history[penalty_offset_];
|
||||
double pij = *penalty_offset;
|
||||
const double wij = std::max(1.0-pij,0.0);
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", " << gm->contact_type << ", " << *gm->xi[1] << ", " << *gm->xj[1] << std::endl;
|
||||
|
||||
// area related calculations
|
||||
double Ac;
|
||||
(*Yflag_offset == 0.0) ? Ac = M_PI*delta*R : Ac = M_PI*((2.0*delta*R - pow(delta,2.0)) + cA/M_PI);
|
||||
if (Ac < 0.0 ) Ac = 0.0;
|
||||
Atot_sum[i] += wij*(Ac - 2.0*M_PI*R*(deltamax_MDR + delta_BULK));
|
||||
Acon1[i] += wij*Ac;
|
||||
|
||||
// bulk force calculation
|
||||
double F_BULK;
|
||||
(delta_BULK <= 0.0) ? F_BULK = 0.0 : F_BULK = (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac;
|
||||
|
||||
//if (atom->tag[i_true] == 4135 && lmp->update->ntimestep > 10935000 && gm->contact_type == 2) {
|
||||
// std::cout << "CS: " << contactSide << ", i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
|
||||
//}
|
||||
|
||||
//if (i == 35 && lmp->update->ntimestep % 1000 == 0) {
|
||||
//double **x = atom->x;
|
||||
//const double xi = x[i_true][0];
|
||||
// std::cout << "i: " << i << ", j: " << j << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", j_true: " << j_true << ", j_tag " << atom->tag[j_true] << ", radi_true: " << radi_true << ", radj_true " << radj_true << ", gm->radi: " << gm->radi << ", gm->radj: " << gm->radj << ", R: " << R << ", Ro: " << Ro << std::endl;
|
||||
//}
|
||||
|
||||
//if (i_true == 35 && lmp->update->ntimestep >= 736002 && lmp->update->ntimestep <= 736005) {
|
||||
// //double **x = atom->x;
|
||||
// //const double xi = x[i_true][0];
|
||||
// std::cout << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", R: " << R << ", Ro: " << Ro << std::endl;
|
||||
//}
|
||||
|
||||
//if ( ((atom->tag[i_true] == 4 && atom->tag[j_true] == 11) || (atom->tag[i_true] == 11 && atom->tag[j_true] == 4)) && lmp->update->ntimestep == 45000) {
|
||||
// std::cout << "CS: " << contactSide << ", contact_type: " << gm->contact_type << ", pair: " << PAIR << ", wall: " << WALL << ", WALLREGION " << WALLREGION << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
|
||||
//}
|
||||
|
||||
//if ( ((atom->tag[i_true] == 4301) || (atom->tag[j_true] == 4076)) && lmp->update->ntimestep > 1016700) {
|
||||
// std::cout << "i: " << i << ", j: " << j << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", j_true: " << j_true << ", j_tag " << atom->tag[j_true] << ", radi_true: " << radi_true << ", radj_true " << radj_true << ", gm->radi: " << gm->radi << ", gm->radj: " << gm->radj << ", R: " << R << ", Ro: " << Ro << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
|
||||
//}
|
||||
|
||||
//if ( ((atom->tag[i_true] == 4) || (atom->tag[j_true] == 4)) && lmp->update->ntimestep == 45000) {
|
||||
// int rank = 0;
|
||||
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
// double **x = atom->x;
|
||||
// const double xi = x[i][0];
|
||||
// const double xj = x[j][0];
|
||||
// const double yi = x[i][1];
|
||||
// const double yj = x[j][1];
|
||||
// const double zi = x[i][2];
|
||||
// const double zj = x[j][2];
|
||||
// std::cout << "CS: " << contactSide << ", rank, " << rank << ", CT: " << gm->contact_type << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", i: " << i << ", j: " << j << ", nlocal: " << atom->nlocal << ", nghost: " << atom->nghost << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << ", R: " << R << ", xi: " << xi << ", xj: " << xj << ", yi: " << yi << ", yj: " << yj << ", zi: " << zi << ", zj: " << zj << std::endl;
|
||||
//}
|
||||
|
||||
//int rank = 0;
|
||||
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
//double **x = atom->x;
|
||||
//const double xi = x[i][0];
|
||||
//const double xj = x[j][0];
|
||||
//const double yi = x[i][1];
|
||||
//const double yj = x[j][1];
|
||||
//const double zi = x[i][2];
|
||||
//const double zj = x[j][2];
|
||||
//const double dis = sqrt(pow((xi-xj),2.0) + pow((yi-yj),2.0) + pow((zi-zj),2.0));
|
||||
//const double delta_test = gm->radi + gm->radj - dis;
|
||||
//if (delta_test < 0.0 && gm->contact_type != 2) {
|
||||
// std::cout << "Particles are not touching but a force is evaluated, CS: " << contactSide << ", rank, " << rank << ", contact_type: " << gm->contact_type << ", pair: " << PAIR << ", wall: " << WALL << ", WALLREGION " << WALLREGION << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", i: " << i << ", j: " << j << ", nlocal: " << atom->nlocal << ", nghost: " << atom->nghost << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << ", R: " << R << ", xi: " << xi << ", xj: " << xj << ", yi: " << yi << ", yj: " << yj << ", zi: " << zi << ", zj: " << zj << std::endl;
|
||||
// std::exit(1);
|
||||
//}
|
||||
|
||||
//if (F_BULK > 0.0) {
|
||||
// std::cout << "F_BULK is: " << F_BULK << std::endl;
|
||||
//}
|
||||
|
||||
//std::cout << delta_BULK << ", " << F_BULK << ", " << (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac << std::endl;
|
||||
|
||||
//int rank = 0;
|
||||
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
//std::cout << "Step: " << lmp->update->ntimestep << ", CS: " << contactSide << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", rank: " << rank << ", gm->delta: " << gm->delta << ", delta: " << delta << ", ddelta: " << ddelta << ", delta_bar: " << ddelta_bar[i] << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << std::endl;
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", " << Vgeo[i] << ", " << Acon0[i] << ", " << Acon1[i] << ", " << Ac << ", " << kappa << " || " << psi[i] << ", " << ddelta_bar[i] << ", " << ddelta << ", " << ddelta_MDR << ", " << ddelta_BULK << ", " << delta << ", " << delta_MDR << ", " << delta_BULK << ", " << F_MDR << ", " << F_BULK << ", " << R << " || " << deltae1D << ", " << A << ", " << B << std::endl;
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", " << (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac << std::endl;
|
||||
|
||||
// total force calculation
|
||||
(contactSide == 0) ? F0 = F_MDR + F_BULK : F1 = F_MDR + F_BULK;
|
||||
|
||||
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << " | " << deltao << ", " << ddelta_bar[i] << ", " << R << ", " << psi[i] << ", " << psi_b << ", " << Ac << " | " << pij << ", " << wij << std::endl;
|
||||
|
||||
// mean surface dipslacement calculation
|
||||
*Ac_offset = wij*Ac;
|
||||
|
||||
// radius update scheme quantity calculation
|
||||
Vcaps[i] += (M_PI/3.0)*pow(delta,2.0)*(3.0*R - delta);
|
||||
|
||||
const double Fntmp = wij*(F_MDR + F_BULK);
|
||||
const double fx = Fntmp*gm->nx[0];
|
||||
const double fy = Fntmp*gm->nx[1];
|
||||
const double fz = Fntmp*gm->nx[2];
|
||||
const double bx = -(Ro - deltao)*gm->nx[0];
|
||||
const double by = -(Ro - deltao)*gm->nx[1];
|
||||
const double bz = -(Ro - deltao)*gm->nx[2];
|
||||
const double eps_bar_contact = (1.0/(3*kappa*Velas[i]))*(fx*bx + fy*by + fz*bz);
|
||||
eps_bar[i] += eps_bar_contact;
|
||||
|
||||
//double **x = atom->x;
|
||||
//const double xi = x[i_true][0];
|
||||
//const double xj = x[j_true][0];
|
||||
//std::cout << i_true << ", " << j_true << ", " << xi << ", " << xj << ", " << gm->nx[0] << ", " << gm->nx[1] << ", " << gm->nx[2] << std::endl;
|
||||
|
||||
//if ( (i == 0 && j == 2 && gm->contact_type == 0) || (i == 2 && j == 0 && gm->contact_type == 0)) {
|
||||
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << delta << ", " << *delta_offset << ", " << (uintptr_t)(delta_offset) << " || " << deltao << ", " << *deltao_offset << ", " << (uintptr_t)(deltao_offset) << " || " << delta_MDR << ", " << *delta_MDR_offset << ", " << (uintptr_t)(delta_MDR_offset) << " || " << *Yflag_offset << ", " << (uintptr_t)(Yflag_offset) << " || " << R << std::endl;
|
||||
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << fx << ", " << fy << ", " << fz << " || " << bx << ", " << by << ", " << bz << ", " << Velas[i] << std::endl;
|
||||
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << eps_bar_contact << ", " << *eps_bar_offset << ", " << (uintptr_t)(eps_bar_offset) << " || " << wij << ", " << ddeltao << ", " << deltao << ", " << delta << ", " << *delta_offset << " || " << Ro << ", " << R << std::endl;
|
||||
//}
|
||||
|
||||
//if () {
|
||||
// std::cout << j << ", " << -Vo*(eps_bar_contact - *eps_bar_offset) - wij*M_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) ) << ", " << -Vo*(eps_bar_contact - *eps_bar_offset) << ", " << wij*M_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) ) << std::endl;
|
||||
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << eps_bar_contact << ", " << *eps_bar_offset << ", " << (uintptr_t)(eps_bar_offset) << " || " << wij << ", " << ddeltao << ", " << deltao << " || " << Ro << ", " << R << std::endl;
|
||||
//}
|
||||
|
||||
|
||||
double desp_bar_contact = eps_bar_contact - *eps_bar_offset; // && desp_bar_contact < 0.0
|
||||
if(delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){
|
||||
const double Vo = (4.0/3.0)*M_PI*pow(Ro,3.0);
|
||||
dRnumerator[i] += -Vo*(eps_bar_contact - *eps_bar_offset) - wij*M_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) );
|
||||
dRdenominator[i] += wij*2.0*M_PI*R*(deltao + R - Ro);
|
||||
|
||||
//if ( (atom->tag[i] == 9) ) {
|
||||
// std::cout << "CT: " << gm->contact_type << ", " << PAIR << "i_tag: " << atom->tag[i] << ", j_tag: " << atom->tag[j] << ", deltae1D: " << deltae1D << ", R: " << R << ", Ro: " << Ro << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", deltao: " << deltao << ", ddeltao: " << ddeltao << ", desp_bar: " << eps_bar_contact - *eps_bar_offset << std::endl;
|
||||
//}
|
||||
}
|
||||
*eps_bar_offset = eps_bar_contact;
|
||||
|
||||
sigmaxx[i] += (1.0/Velas[i])*(fx*bx);
|
||||
sigmayy[i] += (1.0/Velas[i])*(fy*by);
|
||||
sigmazz[i] += (1.0/Velas[i])*(fz*bz);
|
||||
//std::cout << psi_b << ", " << psi[i] << ", " << A << ", " << B << ", " << pY << ", " << amax << " || " << deltao << ", " << delta << ", " << ddelta << ", " << *delta_offset << ", " << ddelta_bar[i] << " || " << delta_MDR << ", " << ddelta_MDR << ", " << *delta_MDR_offset << ", " << deltamax_MDR << " || " << delta_BULK << ", " << ddelta_BULK << ", " << *delta_BULK_offset << " || " << R << " || " << Ac << ", " << *Ac_offset << ", " << Acon0[i] << ", " << Acon1[i] << " || " << F_MDR << ", " << F_BULK << ", " << Vgeo[i] << std::endl;
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", " << gm->radi << ", " << gm->radj << " | " << delta << ", " << F_MDR << " | " << deltae1D << ", " << A << ", " << B << std::endl;
|
||||
|
||||
//if (atom->tag[i] == 1 && lmp->update->ntimestep == 45000) {
|
||||
// double nx;
|
||||
// double ny;
|
||||
// double nz;
|
||||
// if (i == j_true) {
|
||||
// nx = gm->nx[0];
|
||||
// ny = gm->nx[1];
|
||||
// nz = gm->nx[2];
|
||||
// } else {
|
||||
// nx = -gm->nx[0];
|
||||
// ny = -gm->nx[1];
|
||||
// nz = -gm->nx[2];
|
||||
// }
|
||||
// double deltae = deltamax_MDR - *deltap_offset;
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/MPFEM/deformed_particle_shape_data.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
|
||||
// rowDataStream << nx << ", " << ny << ", " << nz << ", " << Ro << ", " << R << ", " << delta << ", " << deltae << ", " << A << ", " << B << ", " << a_na;
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
}
|
||||
|
||||
gm->i = i_true;
|
||||
gm->j = j_true;
|
||||
gm->radi = radi_true;
|
||||
gm->radj = radj_true;
|
||||
|
||||
double * penalty_offset = & history[penalty_offset_];
|
||||
const double pij = *penalty_offset;
|
||||
const double wij = std::max(1.0-pij,0.0);
|
||||
*penalty_offset = 0.0;
|
||||
|
||||
//std::cout << gm->i << ", " << gm->j << ", " << xi << ", " << xj << std::endl;
|
||||
|
||||
// force magnifiers to prevent over penetration
|
||||
double * deltao_offset = & history[deltao_offset_0];
|
||||
//const double wallForceMagnifer = std::exp(10.0*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0;
|
||||
const double wallForceMagnifer = 1.0;
|
||||
|
||||
// assign final force
|
||||
//(gm->contact_type != PAIR) ? F = wij*F0*wallForceMagnifer : F = wij*(F0 + F1)/2; // F = 2*wij*pow(1/F0 + 1/F1,-1);
|
||||
|
||||
if (gm->contact_type != PAIR) {
|
||||
F = wij*F0*wallForceMagnifer;
|
||||
} else {
|
||||
F = wij*(F0 + F1)/2.0;
|
||||
}
|
||||
|
||||
//std::cout << F << ", " << F0 << ", " << F1 << " | " << R0 << ", " << R1 << std::endl;
|
||||
|
||||
// calculate damping force
|
||||
if (F > 0.0) {
|
||||
double Eeff;
|
||||
double Reff;
|
||||
if (gm->contact_type == PAIR) {
|
||||
Eeff = E/(2.0*(1.0-pow(nu,2.0)));
|
||||
Reff = pow((1/gm->radi + 1/gm->radj),-1);
|
||||
} else {
|
||||
Eeff = E/(1.0-pow(nu,2.0));
|
||||
Reff = gm->radi;
|
||||
}
|
||||
const double kn = Eeff*Reff;
|
||||
const double beta = -log(CoR)/sqrt(pow(log(CoR),2.0) + M_PI*M_PI);
|
||||
const double damp_prefactor = beta*sqrt(gm->meff*kn);
|
||||
const double F_DAMP = -damp_prefactor*(gm->vnnr);
|
||||
|
||||
//std:: cout << gm->contact_type << ", " << Eeff << " , " << Reff << ", " << gm->radi << ", " << gm->radj << " || " << kn << ", " << beta << ", " << gm->meff << " || " << F_DAMP << ", " << F << std::endl;
|
||||
F += wij*F_DAMP;
|
||||
}
|
||||
|
||||
//double **x = atom->x;
|
||||
//const double xi = x[gm->i][0];
|
||||
//const double xj = x[gm->j][0];
|
||||
//const double del = 20.0 - abs(xi-xj);
|
||||
|
||||
//if (i_true == 146 && j_true == 152) {
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/rigid_flat_output.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
|
||||
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << dde0 << ", " << dde1 << ", " << (dde0 + dde1) << ", " << F0old << ", " << F1old << ", " << a0 << ", " << a1 << ", " << h0 << ", " << h1 << ", " << k_BULK0 << ", " << k_BULK1 << ", " << h_BULK0 << ", " << h_BULK1;
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
//if (i_true == 0 && j_true == 1) {
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/twoParticleDifferingRadii/rigid_flat_output.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
|
||||
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << dde0 << ", " << dde1 << ", " << (dde0 + dde1) << ", " << F0old << ", " << F1old << ", " << a0 << ", " << a1 << ", " << h0 << ", " << h1 << ", " << k_BULK0 << ", " << k_BULK1 << ", " << h_BULK0 << ", " << h_BULK1 << ", " << R0 << ", " << R1;
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
//if (i_true == 0 && j_true == 1) {
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/twoParticleDifferingRadii/rigid_flat_output_fit.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
|
||||
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << delta01;
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
//if (gm->i == 0 && gm->j == 2) {
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps/sims/compressionSleeve/pairContactsBotCen.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
|
||||
// rowDataStream << del << ", " << F;
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
//std:: cout << "The force F is: " << F << std::endl;
|
||||
|
||||
return F;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void GranSubModNormalMDR::set_fncrit()
|
||||
{
|
||||
Fncrit = fabs(F);
|
||||
}
|
||||
|
||||
//std::cout << sidata.i << ", " << sidata.j << ", " << R << ", " << deltan << ", " << deltao << ", " << dRsums_i[0] << ", " << dRsums_i[1] << ", " << numQuant << std::endl;
|
||||
//std::cout << gm->i << ", " << gm->j << " | " << gm->nx[0] << ", " << gm->nx[1] << ", " << gm->nx[2] << std::endl;
|
||||
//std::cout << "F is: " << F << std::endl;
|
||||
//std::cout << "Ro from particle history is: " << Ro[gm->i] << std::endl;
|
||||
//std::cout << "MDR contact model has been entered." << std::endl;
|
||||
// std::cout << "F is: " << F << std::endl;
|
||||
//std::cout << "gamma > 0.0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
|
||||
//std::cout << "deltae1D <= 0.0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
|
||||
//std::cout << "F_MDR should be > 0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
|
||||
//std::cout << gm->i << ", " << gm->j << " || " << delta << ", " << delta_MDR << ", " << deltamax_MDR << ", " << deltae1D << " || " << A << ", " << B << std::endl;
|
||||
|
||||
//std::cout << i_true << ", " << j_true << std::endl;
|
||||
|
||||
//std::cout << history_index << ", " << history[0] << ", " << history[1] << ", " << history[2] << std::endl;
|
||||
|
||||
// initialize all history variables
|
||||
//double delta_offset;
|
||||
//double deltao_offset;
|
||||
//double delta_MDR_offset;
|
||||
//double delta_BULK_offset;
|
||||
//double deltamax_MDR_offset;
|
||||
//double Yflag;
|
||||
//double deltaY_offset;
|
||||
//double Ac_offset;
|
||||
//double aAdh_offset;
|
||||
//double deltap_offset;
|
||||
//double cA_offset;
|
||||
//double eps_bar_offset;
|
||||
//double wall_contact_flag_offset;
|
||||
@ -19,6 +19,7 @@ GranSubModStyle(hertz,GranSubModNormalHertz,NORMAL);
|
||||
GranSubModStyle(hertz/material,GranSubModNormalHertzMaterial,NORMAL);
|
||||
GranSubModStyle(dmt,GranSubModNormalDMT,NORMAL);
|
||||
GranSubModStyle(jkr,GranSubModNormalJKR,NORMAL);
|
||||
GranSubModStyle(mdr,GranSubModNormalMDR,NORMAL);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
@ -133,6 +134,20 @@ namespace Granular_NS {
|
||||
int mixed_coefficients;
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
class GranSubModNormalMDR : public GranSubModNormal {
|
||||
public:
|
||||
GranSubModNormalMDR(class GranularModel *, class LAMMPS *);
|
||||
void coeffs_to_local() override;
|
||||
double calculate_forces() override;
|
||||
void set_fncrit() override;
|
||||
double psi_b;
|
||||
|
||||
protected:
|
||||
double E, nu, Y, gamma, CoR, F;
|
||||
};
|
||||
|
||||
} // namespace Granular_NS
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
|
||||
@ -17,6 +17,7 @@
|
||||
#include "gran_sub_mod_normal.h"
|
||||
#include "granular_model.h"
|
||||
#include "math_extra.h"
|
||||
#include <iostream>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -46,7 +47,7 @@ GranSubModRollingNone::GranSubModRollingNone(GranularModel *gm, LAMMPS *lmp) :
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
GranSubModRollingSDS::GranSubModRollingSDS(GranularModel *gm, LAMMPS *lmp) :
|
||||
GranSubModRolling(gm, lmp)
|
||||
GranSubModRolling(gm, lmp), k{0.0}, mu{0.0}, gamma{0.0}
|
||||
{
|
||||
num_coeffs = 3;
|
||||
size_history = 3;
|
||||
@ -81,6 +82,8 @@ void GranSubModRollingSDS::calculate_forces()
|
||||
hist_temp[1] = gm->history[rhist1];
|
||||
hist_temp[2] = gm->history[rhist2];
|
||||
|
||||
//std::cout << "Frcrit rolling is: " << Frcrit << std::endl;
|
||||
|
||||
if (gm->history_update) {
|
||||
rolldotn = dot3(hist_temp, gm->nx);
|
||||
|
||||
|
||||
@ -243,7 +243,7 @@ void GranularModel::init()
|
||||
|
||||
// Must have valid normal, damping, and tangential models
|
||||
if (normal_model->name == "none") error->all(FLERR, "Must specify normal granular model");
|
||||
if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model");
|
||||
//if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); TURNED OFF FOR MDR
|
||||
if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model");
|
||||
|
||||
// Twisting, rolling, and heat are optional
|
||||
|
||||
@ -74,6 +74,9 @@ class GranularModel : protected Pointers {
|
||||
int beyond_contact, limit_damping, history_update;
|
||||
ContactType contact_type;
|
||||
|
||||
// Particle indices
|
||||
int i, j, itype, jtype;
|
||||
|
||||
// History variables
|
||||
int size_history, nondefault_history_transfer;
|
||||
double *transfer_history_factor;
|
||||
|
||||
@ -32,11 +32,14 @@
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "group.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "update.h"
|
||||
#include "gran_sub_mod_normal.h"
|
||||
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
@ -78,6 +81,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
|
||||
|
||||
fix_history = nullptr;
|
||||
fix_dummy = dynamic_cast<FixDummy *>(modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY"));
|
||||
|
||||
fix_flag = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -194,11 +199,28 @@ void PairGranular::compute(int eflag, int vflag)
|
||||
jtype = type[j];
|
||||
model = models_list[types_indices[itype][jtype]];
|
||||
|
||||
//std::cout << "Pair Granular: " << i << ", " << j << std::endl;
|
||||
|
||||
// Reset model and copy initial geometric data
|
||||
model->xi = x[i];
|
||||
model->xj = x[j];
|
||||
model->radi = radius[i];
|
||||
model->radj = radius[j];
|
||||
model->i = i;
|
||||
model->j = j;
|
||||
model->itype = itype;
|
||||
model->jtype = jtype;
|
||||
|
||||
/*
|
||||
if (radius[i] < 0.00065 || radius[j] < 0.00065) {
|
||||
std::cout << "Pair Granular: " << i << ", " << j << " || " << radius[i] << ", " << radius[j] << std::endl;
|
||||
}
|
||||
|
||||
if (std::isnan(radius[i]) || std::isnan(radius[j])) {
|
||||
std::cout << "Pair Granular: " << i << ", " << j << " || " << radius[i] << ", " << radius[j] << std::endl;
|
||||
}
|
||||
*/
|
||||
|
||||
if (use_history) model->touch = touch[jj];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
@ -468,6 +490,67 @@ void PairGranular::init_style()
|
||||
if (!fix_history) error->all(FLERR,"Could not find pair fix neigh history ID");
|
||||
}
|
||||
|
||||
/*
|
||||
// Example
|
||||
//Store persistent per atom quantities
|
||||
if (! fix_flag) {
|
||||
int tmp1, tmp2;
|
||||
id_fix = "BOND_BPM_PLASTIC_FIX_PROP_ATOM";
|
||||
modify->add_fix(
|
||||
fmt::format("{} all property/atom d_plastic d_plastic_temp d_viscous_temp d_plastic_heat d_viscous_heat ghost yes", id_fix));
|
||||
index_plastic = atom->find_custom("plastic",tmp1,tmp2);
|
||||
index_pq = atom->find_custom("plastic_heat",tmp1,tmp2);
|
||||
index_pt = atom->find_custom("plastic_temp",tmp1,tmp2);
|
||||
index_vq = atom->find_custom("viscous_heat",tmp1,tmp2);
|
||||
index_vt = atom->find_custom("viscous_temp",tmp1,tmp2);
|
||||
fix_flag = 1;
|
||||
}
|
||||
*/
|
||||
|
||||
if (model->normal_model->name == "mdr") {
|
||||
|
||||
std::cout << "MDR history variables have been initialized" << std::endl;
|
||||
|
||||
// FOR MDR CONTACT MODEL
|
||||
//Store persistent per atom quantities
|
||||
if (! fix_flag) {
|
||||
int tmp1, tmp2;
|
||||
const char * id_fix = "MDR_PARTICLE_HISTORY_VARIABLES";
|
||||
modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_psi_b d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_contacts d_adhesive_length ghost yes", id_fix));
|
||||
// d2_volSums 4 --> allows an array of 4 to defined.
|
||||
index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
|
||||
index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
|
||||
index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
|
||||
index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
|
||||
index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
|
||||
index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
|
||||
index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
|
||||
index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
|
||||
index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
|
||||
index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
|
||||
index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
|
||||
index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
|
||||
index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
|
||||
index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THE INPUT SCRIPT
|
||||
index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); // flag to check if history variables have been initialized
|
||||
index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation
|
||||
index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
|
||||
index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
|
||||
index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
|
||||
index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
|
||||
|
||||
std::cout << "MDR history variables have been initialized 2" << ", " << index_Ro << std::endl;
|
||||
|
||||
//index_volSums = atom->find_custom("volSums",tmp1,tmp2);
|
||||
|
||||
// Initiate MDR radius update fix
|
||||
modify->add_fix("fix_mdr_radius_update all mdr/radius/update");
|
||||
modify->add_fix("fix_mdr_mean_surf_disp all mdr/mean/surf/disp");
|
||||
|
||||
fix_flag = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// check for FixFreeze and set freeze_group_bit
|
||||
|
||||
auto fixlist = modify->get_fix_by_style("^freeze");
|
||||
@ -703,6 +786,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
||||
|
||||
class GranularModel* model = models_list[types_indices[itype][jtype]];
|
||||
|
||||
|
||||
|
||||
// Reset model and copy initial geometric data
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
@ -925,3 +1010,8 @@ void PairGranular::prune_models()
|
||||
nmodels -= 1;
|
||||
}
|
||||
}
|
||||
|
||||
size_t PairGranular::get_size_history() const
|
||||
{
|
||||
return size_history;
|
||||
}
|
||||
|
||||
@ -46,6 +46,16 @@ class PairGranular : public Pair {
|
||||
double memory_usage() override;
|
||||
double atom2cut(int) override;
|
||||
double radii2cut(double, double) override;
|
||||
size_t get_size_history() const;
|
||||
|
||||
// granular models
|
||||
// MOVED HERE FROM PRIVATE FOR MDR MODEL
|
||||
class Granular_NS::GranularModel** models_list;
|
||||
int **types_indices;
|
||||
int nmodels, maxmodels;
|
||||
|
||||
class FixStoreAtom * fix_store;
|
||||
|
||||
|
||||
protected:
|
||||
int freeze_group_bit;
|
||||
@ -59,6 +69,29 @@ class PairGranular : public Pair {
|
||||
class FixDummy *fix_dummy;
|
||||
class FixNeighHistory *fix_history;
|
||||
|
||||
// MDR particle history variables
|
||||
int fix_flag;
|
||||
int index_Ro;
|
||||
int index_Vcaps;
|
||||
int index_Vgeo;
|
||||
int index_Velas;
|
||||
int index_eps_bar;
|
||||
int index_dRnumerator;
|
||||
int index_dRdenominator;
|
||||
int index_Acon0;
|
||||
int index_Acon1;
|
||||
int index_Atot;
|
||||
int index_Atot_sum;
|
||||
int index_ddelta_bar;
|
||||
int index_psi;
|
||||
int index_psi_b;
|
||||
int index_history_setup_flag;
|
||||
int index_sigmaxx;
|
||||
int index_sigmayy;
|
||||
int index_sigmazz;
|
||||
int index_contacts;
|
||||
int index_adhesive_length;
|
||||
|
||||
// storage of rigid body masses for use in granular interactions
|
||||
|
||||
class Fix *fix_rigid; // ptr to rigid body fix, null pointer if none
|
||||
@ -73,11 +106,6 @@ class PairGranular : public Pair {
|
||||
int size_history;
|
||||
int heat_flag;
|
||||
|
||||
// granular models
|
||||
int nmodels, maxmodels;
|
||||
class Granular_NS::GranularModel **models_list;
|
||||
int **types_indices;
|
||||
|
||||
// optional user-specified global cutoff, per-type user-specified cutoffs
|
||||
double **cutoff_type;
|
||||
double cutoff_global;
|
||||
|
||||
@ -40,6 +40,7 @@
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
@ -54,7 +55,7 @@ enum { SHIFT_NO, SHIFT_YES, SHIFT_POSSIBLE };
|
||||
|
||||
static constexpr double EINERTIA = 0.2; // moment of inertia prefactor for ellipsoid
|
||||
|
||||
static constexpr int ATOMPERBIN = 30;
|
||||
static constexpr int ATOMPERBIN = 100;
|
||||
static constexpr double BIG = 1.0e20;
|
||||
static constexpr int VBINSIZE = 5;
|
||||
static constexpr double TOLERANCE = 0.00001;
|
||||
@ -1321,6 +1322,7 @@ void FixSRD::collisions_single()
|
||||
#endif
|
||||
|
||||
if (t_remain > dt) {
|
||||
print_collision(i, j, ibounce, t_remain, dt, xscoll, xbcoll, norm, type);
|
||||
ninside++;
|
||||
if (insideflag == INSIDE_ERROR || insideflag == INSIDE_WARN) {
|
||||
std::string mesg;
|
||||
@ -2603,14 +2605,17 @@ void FixSRD::parameterize()
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) {
|
||||
if (mass_srd == 0.0)
|
||||
if (mass_srd == 0.0) {
|
||||
mass_srd = rmass[i];
|
||||
else if (rmass[i] != mass_srd)
|
||||
std::cout << "mass if #1: " << rmass[i] << std::endl;
|
||||
std::cout << "srd radius is: " << radius[i] << std::endl;
|
||||
} else if (rmass[i] != mass_srd)
|
||||
flag = 1;
|
||||
} else {
|
||||
if (mass_srd == 0.0)
|
||||
if (mass_srd == 0.0) {
|
||||
mass_srd = mass[type[i]];
|
||||
else if (mass[type[i]] != mass_srd)
|
||||
std::cout << "mass if #2: " << mass[type[i]] << std::endl;
|
||||
} else if (mass[type[i]] != mass_srd)
|
||||
flag = 1;
|
||||
}
|
||||
}
|
||||
@ -2712,6 +2717,7 @@ void FixSRD::parameterize()
|
||||
if (dimension == 3) {
|
||||
volsrd = (domain->xprd * domain->yprd * domain->zprd) - volbig;
|
||||
density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd * domain->zprd - volbig);
|
||||
std::cout << nsrd << ", " << mass_srd << ", " << domain->xprd << ", " << domain->yprd << ", " << domain->zprd << ", " << volbig << std::endl;
|
||||
} else {
|
||||
volsrd = (domain->xprd * domain->yprd) - volbig;
|
||||
density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd - volbig);
|
||||
@ -2773,7 +2779,7 @@ void FixSRD::parameterize()
|
||||
gridsrd, binsize3x, binsize3y, binsize3z);
|
||||
mesg += fmt::format(" SRD per actual grid cell = {:.8}\n", srd_per_cell);
|
||||
mesg += fmt::format(" SRD viscosity = {:.8}\n", viscosity);
|
||||
mesg += fmt::format(" big/SRD mass density ratio = {:.8}\n", mdratio);
|
||||
mesg += fmt::format(" big/SRD mass density ratio = {:.8} | big density = {:.8}, small density = {:.8}\n", mdratio,density_big,density_srd);
|
||||
utils::logmesg(lmp, mesg);
|
||||
}
|
||||
|
||||
|
||||
25
src/csv_writer.h
Normal file
25
src/csv_writer.h
Normal file
@ -0,0 +1,25 @@
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
|
||||
class CSVWriter {
|
||||
public:
|
||||
CSVWriter(const std::string& filename) : filename_(filename) {}
|
||||
|
||||
void writeRow(const std::string& data) {
|
||||
std::ofstream file;
|
||||
// Use the append mode to add data to the end of the file if it exists
|
||||
file.open(filename_, std::ios::out | std::ios::app);
|
||||
|
||||
if (!file.is_open()) {
|
||||
std::cerr << "Failed to open file: " << filename_ << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
file << data << std::endl;
|
||||
file.close();
|
||||
}
|
||||
|
||||
private:
|
||||
std::string filename_;
|
||||
};
|
||||
@ -431,6 +431,8 @@ FixAveChunk::~FixAveChunk()
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
// FOR MDR
|
||||
|
||||
int FixAveChunk::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
@ -499,6 +501,8 @@ void FixAveChunk::setup(int /*vflag*/)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
// FOR MDR, DO WHATEVER YOUR FIX NEEDS TO DO.
|
||||
|
||||
void FixAveChunk::end_of_step()
|
||||
{
|
||||
int i,j,m,index;
|
||||
|
||||
@ -31,7 +31,7 @@ class FixAveChunk : public Fix {
|
||||
int setmask() override;
|
||||
void init() override;
|
||||
void setup(int) override;
|
||||
void end_of_step() override;
|
||||
void end_of_step() override; // FOR MDR
|
||||
double compute_array(int, int) override;
|
||||
double memory_usage() override;
|
||||
|
||||
|
||||
539
src/fix_mdr_mean_surf_disp.cpp
Normal file
539
src/fix_mdr_mean_surf_disp.cpp
Normal file
@ -0,0 +1,539 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors:
|
||||
William Zunker (MIT), Sachith Dunatunga (MIT),
|
||||
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
|
||||
----------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_mdr_mean_surf_disp.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "input.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "variable.h"
|
||||
#include "fix_neigh_history.h"
|
||||
#include "pair.h"
|
||||
#include "pair_granular.h"
|
||||
#include "granular_model.h"
|
||||
#include "neigh_list.h"
|
||||
#include "region.h"
|
||||
#include "update.h"
|
||||
#include "fix_wall_gran_region.h"
|
||||
#include "comm.h"
|
||||
#include <iostream>
|
||||
#include <algorithm>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace Granular_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixMDRmeanSurfDisp::FixMDRmeanSurfDisp(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
comm_forward = 1; // value needs to match number of values you communicate
|
||||
}
|
||||
|
||||
// FOR MDR
|
||||
|
||||
int FixMDRmeanSurfDisp::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
void FixMDRmeanSurfDisp::setup_pre_force(int /*vflag*/)
|
||||
{
|
||||
int tmp1, tmp2;
|
||||
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
|
||||
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
|
||||
Acon0 = atom->dvector[index_Acon0];
|
||||
ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
|
||||
pre_force(0);
|
||||
}
|
||||
|
||||
int FixMDRmeanSurfDisp::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
|
||||
{
|
||||
int m = 0;
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
buf[m++] = ddelta_bar[j];
|
||||
//buf[m++] = Acon0[j];
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
void FixMDRmeanSurfDisp::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int m = 0;
|
||||
int last = first + n;
|
||||
for (int i = first; i < last; i++) {
|
||||
ddelta_bar[i] = buf[m++];
|
||||
//Acon0[i] = buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
void FixMDRmeanSurfDisp::pre_force(int)
|
||||
{
|
||||
FixNeighHistory * fix_history = dynamic_cast<FixNeighHistory *>(modify->get_fix_by_id("NEIGH_HISTORY_GRANULAR"));
|
||||
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
|
||||
NeighList * list = pair->list;
|
||||
|
||||
const int size_history = pair->get_size_history();
|
||||
|
||||
//std::cout << " " << std::endl;
|
||||
//std::cout << "New Step" << std::endl;
|
||||
|
||||
{
|
||||
int i,j,k,lv1,ii,jj,inum,jnum,itype,jtype,ktype;
|
||||
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int *touch,**firsttouch;
|
||||
double *history_ij,*history_ik,*history_jk,*history_kj,*allhistory,*allhistory_j,*allhistory_k,**firsthistory;
|
||||
|
||||
bool touchflag = false;
|
||||
|
||||
//class GranularModel* model;
|
||||
//class GranularModel** models_list = pair->models_list;
|
||||
//int ** types_indices = pair->types_indices;
|
||||
|
||||
double **x = atom->x;
|
||||
int *type = atom->type;
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
firsttouch = fix_history->firstflag;
|
||||
firsthistory = fix_history->firstvalue;
|
||||
|
||||
// contact penalty calculation
|
||||
for (int ii = 0; ii < inum; ii++) {
|
||||
const int i = ilist[ii];
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
allhistory = firsthistory[i];
|
||||
double radi = radius[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
jtype = type[j];
|
||||
double radj = radius[j];
|
||||
const double delx_ij = x[j][0] - xtmp;
|
||||
const double dely_ij = x[j][1] - ytmp;
|
||||
const double delz_ij = x[j][2] - ztmp;
|
||||
const double rsq_ij = delx_ij*delx_ij + dely_ij*dely_ij + delz_ij*delz_ij;
|
||||
const double r_ij = sqrt(rsq_ij);
|
||||
const double rinv_ij = 1.0/r_ij;
|
||||
const double radsum_ij = radi + radj;
|
||||
const double deltan_ij = radsum_ij - r_ij;
|
||||
if (deltan_ij >= 0.0) {
|
||||
for (int kk = jj; kk < jnum; kk++) {
|
||||
k = jlist[kk];
|
||||
k &= NEIGHMASK;
|
||||
ktype = type[k];
|
||||
if (kk != jj) {
|
||||
const double delx_ik = x[k][0] - xtmp;
|
||||
const double dely_ik = x[k][1] - ytmp;
|
||||
const double delz_ik = x[k][2] - ztmp;
|
||||
const double rsq_ik = delx_ik*delx_ik + dely_ik*dely_ik + delz_ik*delz_ik;
|
||||
const double r_ik = sqrt(rsq_ik);
|
||||
const double rinv_ik = 1.0/r_ik;
|
||||
const double radk = radius[k];
|
||||
const double radsum_ik = radi + radk;
|
||||
const double deltan_ik = radsum_ik - r_ik;
|
||||
const double delx_jk = x[k][0] - x[j][0];
|
||||
const double dely_jk = x[k][1] - x[j][1];
|
||||
const double delz_jk = x[k][2] - x[j][2];
|
||||
const double rsq_jk = delx_jk*delx_jk + dely_jk*dely_jk + delz_jk*delz_jk;
|
||||
const double r_jk = sqrt(rsq_jk);
|
||||
const double rinv_jk = 1.0/r_jk;
|
||||
const double radsum_jk = radj + radk;
|
||||
const double deltan_jk = radsum_jk - r_jk;
|
||||
if (deltan_ik >= 0.0 && deltan_jk >= 0.0) {
|
||||
|
||||
// pull ij history
|
||||
history_ij = &allhistory[size_history * jj];
|
||||
double * pij = &history_ij[22]; // penalty for contact i and j
|
||||
|
||||
// pull ik history
|
||||
history_ik = &allhistory[size_history * kk];
|
||||
double * pik = &history_ik[22]; // penalty for contact i and k
|
||||
|
||||
// we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa, so we need to manually search to figure out the owner
|
||||
// check if k is in the neighbor list of j
|
||||
double * pjk = NULL;
|
||||
int * const jklist = firstneigh[j];
|
||||
const int jknum = numneigh[j];
|
||||
for (int jk = 0; jk < jknum; jk++) {
|
||||
const int kneigh = jklist[jk] & NEIGHMASK;
|
||||
if (k == kneigh) {
|
||||
allhistory_j = firsthistory[j];
|
||||
history_jk = &allhistory_j[size_history * jk];
|
||||
pjk = &history_jk[22]; // penalty for contact j and k
|
||||
//int rank = 0;
|
||||
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
//std::cout << "Print 183 pjk: " << rank << ", " << pjk << std::endl;
|
||||
//std::cout << "Print 183 pjk[0]: " << rank << ", " << pjk[0] << std::endl;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// check if j is in the neighbor list of k
|
||||
if (pjk == NULL) {
|
||||
int * const kjlist = firstneigh[k];
|
||||
const int kjnum = numneigh[k];
|
||||
for (int kj = 0; kj < kjnum; kj++) {
|
||||
const int jneigh = kjlist[kj] & NEIGHMASK;
|
||||
if (j == jneigh) {
|
||||
allhistory_k = firsthistory[k];
|
||||
history_kj = &allhistory_k[size_history * kj];
|
||||
pjk = &history_kj[22]; // penalty for contact j and k
|
||||
//int rank = 0;
|
||||
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
//std::cout << "Print 198 pjk: " << rank << ", " << pjk << std::endl;
|
||||
//std::cout << "Print 198 pjk[0]: " << rank << ", " << pjk[0] << std::endl;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//std::cout << "Print: " << __LINE__ << std::endl;
|
||||
std::vector<double> distances = {r_ij,r_ik,r_jk};
|
||||
auto maxElement = std::max_element(distances.begin(), distances.end());
|
||||
double maxValue = *maxElement;
|
||||
int maxIndex = std::distance(distances.begin(), maxElement);
|
||||
if (maxIndex == 0) { // the central particle is k
|
||||
const double enx_ki = -delx_ik * rinv_ik;
|
||||
const double eny_ki = -dely_ik * rinv_ik;
|
||||
const double enz_ki = -delz_ik * rinv_ik;
|
||||
const double enx_kj = -delx_jk * rinv_jk;
|
||||
const double eny_kj = -dely_jk * rinv_jk;
|
||||
const double enz_kj = -delz_jk * rinv_jk;
|
||||
const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj);
|
||||
pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
|
||||
} else if (maxIndex == 1) { // the central particle is j
|
||||
const double enx_ji = -delx_ij * rinv_ij;
|
||||
const double eny_ji = -dely_ij * rinv_ij;
|
||||
const double enz_ji = -delz_ij * rinv_ij;
|
||||
const double enx_jk = delx_jk * rinv_jk;
|
||||
const double eny_jk = dely_jk * rinv_jk;
|
||||
const double enz_jk = delz_jk * rinv_jk;
|
||||
const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk);
|
||||
pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
|
||||
} else { // the central particle is i
|
||||
if (j < atom->nlocal || k < atom->nlocal) {
|
||||
const double enx_ij = delx_ij * rinv_ij;
|
||||
const double eny_ij = dely_ij * rinv_ij;
|
||||
const double enz_ij = delz_ij * rinv_ij;
|
||||
const double enx_ik = delx_ik * rinv_ik;
|
||||
const double eny_ik = dely_ik * rinv_ik;
|
||||
const double enz_ik = delz_ik * rinv_ik;
|
||||
const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik);
|
||||
//std::cout << "Print: " << __LINE__ << std::endl;
|
||||
pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
|
||||
//std::cout << "Print: " << __LINE__ << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
{
|
||||
int i,j,k,ii,jj,inum,jnum,itype,jtype;
|
||||
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int *touch,**firsttouch;
|
||||
double *history,*allhistory,**firsthistory;
|
||||
|
||||
bool touchflag = false;
|
||||
|
||||
class GranularModel* model;
|
||||
class GranularModel** models_list = pair->models_list;
|
||||
int ** types_indices = pair->types_indices;
|
||||
|
||||
double **x = atom->x;
|
||||
int *type = atom->type;
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
firsttouch = fix_history->firstflag;
|
||||
firsthistory = fix_history->firstvalue;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itype = type[i];
|
||||
touch = firsttouch[i];
|
||||
allhistory = firsthistory[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
jtype = type[j];
|
||||
model = models_list[types_indices[itype][jtype]];
|
||||
|
||||
// Reset model and copy initial geometric data
|
||||
model->xi = x[i];
|
||||
model->xj = x[j];
|
||||
model->radi = radius[i];
|
||||
model->radj = radius[j];
|
||||
model->i = i;
|
||||
model->j = j;
|
||||
model->touch = touch[jj];
|
||||
touchflag = model->check_contact();
|
||||
|
||||
// is it necessary to clear the history here???
|
||||
if (!touchflag) {
|
||||
touch[jj] = 0;
|
||||
history = &allhistory[size_history * jj];
|
||||
for (k = 0; k < size_history; k++) history[k] = 0.0;
|
||||
continue;
|
||||
}
|
||||
|
||||
touch[jj] = 1;
|
||||
|
||||
history = &allhistory[size_history * jj];
|
||||
model->history = history;
|
||||
|
||||
//const double pij = history[22];
|
||||
//const double wij = std::max(1.0-pij,0.0);
|
||||
const double delta = model->radsum - sqrt(model->rsq);
|
||||
|
||||
const int delta_offset_0 = 0; // apparent overlap
|
||||
const int delta_offset_1 = 1;
|
||||
const int Ac_offset_0 = 18; // contact area
|
||||
const int Ac_offset_1 = 19;
|
||||
const int deltamax_offset_ = 23;
|
||||
const int deltap_offset_0 = 24;
|
||||
const int deltap_offset_1 = 25;
|
||||
|
||||
double deltamax = history[deltamax_offset_];
|
||||
double deltap0 = history[deltap_offset_0];
|
||||
double deltap1 = history[deltap_offset_1];
|
||||
|
||||
if (delta > deltamax) deltamax = delta;
|
||||
|
||||
double delta0old = history[delta_offset_0];
|
||||
double delta1old = history[delta_offset_1];
|
||||
|
||||
int i0;
|
||||
int i1;
|
||||
if (atom->tag[i] > atom->tag[j]) {
|
||||
i0 = i;
|
||||
i1 = j;
|
||||
} else {
|
||||
i0 = j;
|
||||
i1 = i;
|
||||
}
|
||||
|
||||
int rank = 0;
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
int i_ghost;
|
||||
int j_ghost;
|
||||
(i >= atom->nlocal) ? i_ghost = 1 : i_ghost = 0;
|
||||
(j >= atom->nlocal) ? j_ghost = 1 : j_ghost = 0;
|
||||
if (i_ghost == 1) {
|
||||
std::cout << "rank: " << rank << ", i: " << i << ", j: " << j << ", i is ghost? " << i_ghost << ", j is ghost? " << j_ghost << ", nlocal: " << atom->nlocal << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << std::endl;
|
||||
}
|
||||
|
||||
//int i0 = std::max(i,j);
|
||||
//int i1 = std::min(i,j);
|
||||
double R0 = radius[i0];
|
||||
double R1 = radius[i1];
|
||||
|
||||
double delta_geo0;
|
||||
double delta_geo1;
|
||||
double deltaOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1));
|
||||
double deltaOpt2 = deltamax*(deltamax - 2.0*R0)/(2.0*(deltamax - R0 - R1));
|
||||
(R0 < R1) ? delta_geo0 = MAX(deltaOpt1,deltaOpt2) : delta_geo0 = MIN(deltaOpt1,deltaOpt2);
|
||||
(R0 < R1) ? delta_geo1 = MIN(deltaOpt1,deltaOpt2) : delta_geo1 = MAX(deltaOpt1,deltaOpt2);
|
||||
|
||||
double overlap_limit = 0.75;
|
||||
|
||||
if (delta_geo0/R0 > overlap_limit) {
|
||||
delta_geo0 = R0*overlap_limit;
|
||||
delta_geo1 = deltamax - delta_geo0;
|
||||
} else if (delta_geo1/R1 > overlap_limit) {
|
||||
delta_geo1 = R1*overlap_limit;
|
||||
delta_geo0 = deltamax - delta_geo1;
|
||||
}
|
||||
|
||||
double deltap = deltap0 + deltap1;
|
||||
|
||||
double delta0 = delta_geo0 + (deltap0 - delta_geo0)/(deltap - deltamax)*(delta-deltamax);
|
||||
double delta1 = delta_geo1 + (deltap1 - delta_geo1)/(deltap - deltamax)*(delta-deltamax);
|
||||
|
||||
double ddel0 = delta0 - delta0old;
|
||||
double ddel1 = delta1 - delta1old;
|
||||
|
||||
if (Acon0[i0] != 0.0) {
|
||||
const double Ac_offset0 = history[Ac_offset_0];
|
||||
ddelta_bar[i0] += Ac_offset0/Acon0[i0]*ddel0; // Multiply by 0.5 since displacement is shared equally between deformable particles.
|
||||
}
|
||||
|
||||
if (Acon0[i1] != 0.0) {
|
||||
const double Ac_offset1 = history[Ac_offset_1];
|
||||
ddelta_bar[i1] += Ac_offset1/Acon0[i1]*ddel1;
|
||||
}
|
||||
|
||||
//int rank = 0;
|
||||
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
//std::cout << "delta_bar calc: Step: " << lmp->update->ntimestep << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", i: " << i << ", j: " << j << ", rank: " << rank << ", Ac_offset0: " << history[Ac_offset_0] << ", Ac_offset1: " << history[Ac_offset_1] << ", Acon[i0]: " << Acon0[i0] << ", Acon[i1]: " << Acon0[i1] << ", ddel0: " << ddel0 << ", ddel1: " << ddel1 << ", ddelta_bar[i0]: " << ddelta_bar[i0] << ", ddelta_bar[i1]: " << ddelta_bar[i1] << std::endl;
|
||||
|
||||
//if (Acon0[j] != 0.0) {
|
||||
// const double delta_offset0 = history[0];
|
||||
// const double ddelta = delta/2.0 - delta_offset0; // Divide by 2.0 since we are storing 1/2 deltan in main MDR script
|
||||
// const double Ac_offset0 = history[18];
|
||||
// ddelta_bar[j] += Ac_offset0/Acon0[j]*ddelta; // Multiply by 0.5 since displacement is shared equally between deformable particles.
|
||||
//}
|
||||
//
|
||||
//if (Acon0[i] != 0.0) {
|
||||
// const double delta_offset1 = history[1];
|
||||
// const double ddelta = delta/2.0 - delta_offset1; // Divide by 2.0 since we are storing 1/2 deltan in main MDR script
|
||||
// const double Ac_offset1 = history[19];
|
||||
// ddelta_bar[i] += Ac_offset1/Acon0[i]*ddelta;
|
||||
//}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
auto fix_list = modify->get_fix_by_style("wall/gran/region");
|
||||
|
||||
for (int w = 0; w < fix_list.size(); w++) {
|
||||
|
||||
FixWallGranRegion* fix = dynamic_cast<FixWallGranRegion*>(fix_list[w]);
|
||||
GranularModel * model = fix->model;
|
||||
Region * region = fix->region;
|
||||
|
||||
{
|
||||
int i, m, nc, iwall;
|
||||
double vwall[3];
|
||||
bool touchflag = false;
|
||||
|
||||
int history_update = 1;
|
||||
model->history_update = history_update;
|
||||
|
||||
int regiondynamic = region->dynamic_check();
|
||||
if (!regiondynamic) vwall[0] = vwall[1] = vwall[2] = 0.0;
|
||||
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
if (regiondynamic) {
|
||||
region->prematch();
|
||||
region->set_velocity();
|
||||
}
|
||||
|
||||
if (fix->peratom_flag) fix->clear_stored_contacts();
|
||||
|
||||
model->radj = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
if (! region->match(x[i][0], x[i][1], x[i][2])) continue;
|
||||
|
||||
nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i] + model->pulloff_distance(radius[i], 0.0));
|
||||
|
||||
if (nc == 0) {
|
||||
fix->ncontact[i] = 0;
|
||||
continue;
|
||||
}
|
||||
if (nc == 1) {
|
||||
fix->c2r[0] = 0;
|
||||
iwall = region->contact[0].iwall;
|
||||
if (fix->ncontact[i] == 0) {
|
||||
fix->ncontact[i] = 1;
|
||||
fix->walls[i][0] = iwall;
|
||||
for (m = 0; m < size_history; m++) fix->history_many[i][0][m] = 0.0;
|
||||
} else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0])
|
||||
fix->update_contacts(i, nc);
|
||||
} else
|
||||
fix->update_contacts(i, nc);
|
||||
|
||||
|
||||
// process current contacts
|
||||
for (int ic = 0; ic < nc; ic++) {
|
||||
|
||||
// Reset model and copy initial geometric data
|
||||
model->dx[0] = region->contact[ic].delx;
|
||||
model->dx[1] = region->contact[ic].dely;
|
||||
model->dx[2] = region->contact[ic].delz;
|
||||
model->radi = radius[i];
|
||||
model->radj = region->contact[ic].radius;
|
||||
model->r = region->contact[ic].r;
|
||||
|
||||
if (model->beyond_contact) model->touch = fix->history_many[i][fix->c2r[ic]][0];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
|
||||
const double wij = 1.0;
|
||||
|
||||
if (Acon0[i] != 0.0) {
|
||||
const double delta = model->radsum - model->r;
|
||||
const double delta_offset0 = fix->history_many[i][fix->c2r[ic]][0];
|
||||
const double ddelta = delta - delta_offset0;
|
||||
const double Ac_offset0 = fix->history_many[i][fix->c2r[ic]][18];
|
||||
ddelta_bar[i] += wij*Ac_offset0/Acon0[i]*ddelta; // Multiply by 0.5 since displacement is shared equally between deformable particles.
|
||||
//std::cout << delta << ", " << delta_offset0 << " || " << Ac_offset0 << ", " << Acon0[i] << ", " << ddelta << ", " << ddelta_bar[i] << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm->forward_comm(this);
|
||||
|
||||
//and the function delcarations in the header:
|
||||
|
||||
//int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
//void unpack_forward_comm(int, int, double *) override;
|
||||
|
||||
//Then the methods would look like::
|
||||
|
||||
//where comm_stage is a public flag to control hich quantity is being communicated
|
||||
|
||||
}
|
||||
|
||||
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl;
|
||||
//std::cout << "Fix radius update setup has been entered !!!" << std::endl;
|
||||
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;
|
||||
46
src/fix_mdr_mean_surf_disp.h
Normal file
46
src/fix_mdr_mean_surf_disp.h
Normal file
@ -0,0 +1,46 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(mdr/mean/surf/disp,FixMDRmeanSurfDisp);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_MDR_MEAN_SURF_DISP_H
|
||||
#define LMP_FIX_MDR_MEAN_SURF_DISP_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixMDRmeanSurfDisp : public Fix {
|
||||
public:
|
||||
double * Acon0;
|
||||
double * ddelta_bar;
|
||||
|
||||
FixMDRmeanSurfDisp(class LAMMPS *, int, char **);
|
||||
int setmask() override;
|
||||
void setup_pre_force(int) override;
|
||||
void pre_force(int) override; // FOR MDR
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
384
src/fix_mdr_radius_update.cpp
Normal file
384
src/fix_mdr_radius_update.cpp
Normal file
@ -0,0 +1,384 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors:
|
||||
William Zunker (MIT), Sachith Dunatunga (MIT),
|
||||
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
|
||||
----------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_mdr_radius_update.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "input.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "variable.h"
|
||||
#include <iostream>
|
||||
#include "csv_writer.h"
|
||||
#include "granular_model.h"
|
||||
#include "pair_granular.h"
|
||||
#include "pair.h"
|
||||
#include "gran_sub_mod_normal.h"
|
||||
#include "update.h"
|
||||
#include "comm.h"
|
||||
#include <iomanip>
|
||||
#include <sstream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixMDRradiusUpdate::FixMDRradiusUpdate(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
comm_forward = 20; // value needs to match number of values you communicate
|
||||
}
|
||||
|
||||
// FOR MDR
|
||||
|
||||
int FixMDRradiusUpdate::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_FORCE | END_OF_STEP;
|
||||
return mask;
|
||||
}
|
||||
|
||||
void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/)
|
||||
{
|
||||
int tmp1, tmp2;
|
||||
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
|
||||
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
|
||||
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
|
||||
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
|
||||
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
|
||||
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
|
||||
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
|
||||
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
|
||||
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
|
||||
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
|
||||
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
|
||||
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
|
||||
int index_psi = atom->find_custom("psi",tmp1,tmp2);
|
||||
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
|
||||
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
|
||||
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
|
||||
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
|
||||
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
|
||||
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
|
||||
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
|
||||
Ro = atom->dvector[index_Ro];
|
||||
Vgeo = atom->dvector[index_Vgeo];
|
||||
Velas = atom->dvector[index_Velas];
|
||||
Vcaps = atom->dvector[index_Vcaps];
|
||||
eps_bar = atom->dvector[index_eps_bar];
|
||||
dRnumerator = atom->dvector[index_dRnumerator];
|
||||
dRdenominator = atom->dvector[index_dRdenominator];
|
||||
Acon0 = atom->dvector[index_Acon0];
|
||||
Acon1 = atom->dvector[index_Acon1];
|
||||
Atot = atom->dvector[index_Atot];
|
||||
Atot_sum = atom->dvector[index_Atot_sum];
|
||||
ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
psi = atom->dvector[index_psi];
|
||||
psi_b = atom->dvector[index_psi_b];
|
||||
sigmaxx = atom->dvector[index_sigmaxx];
|
||||
sigmayy = atom->dvector[index_sigmayy];
|
||||
sigmazz = atom->dvector[index_sigmazz];
|
||||
history_setup_flag = atom->dvector[index_history_setup_flag];
|
||||
contacts = atom->dvector[index_contacts];
|
||||
adhesive_length = atom->dvector[index_adhesive_length];
|
||||
|
||||
pre_force(0);
|
||||
}
|
||||
|
||||
void FixMDRradiusUpdate::setup(int /*vflag*/)
|
||||
{
|
||||
int tmp1, tmp2;
|
||||
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
|
||||
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
|
||||
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
|
||||
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
|
||||
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
|
||||
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
|
||||
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
|
||||
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
|
||||
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
|
||||
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
|
||||
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
|
||||
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
|
||||
int index_psi = atom->find_custom("psi",tmp1,tmp2);
|
||||
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
|
||||
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
|
||||
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
|
||||
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
|
||||
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
|
||||
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
|
||||
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
|
||||
Ro = atom->dvector[index_Ro];
|
||||
Vgeo = atom->dvector[index_Vgeo];
|
||||
Velas = atom->dvector[index_Velas];
|
||||
Vcaps = atom->dvector[index_Vcaps];
|
||||
eps_bar = atom->dvector[index_eps_bar];
|
||||
dRnumerator = atom->dvector[index_dRnumerator];
|
||||
dRdenominator = atom->dvector[index_dRdenominator];
|
||||
Acon0 = atom->dvector[index_Acon0];
|
||||
Acon1 = atom->dvector[index_Acon1];
|
||||
Atot = atom->dvector[index_Atot];
|
||||
Atot_sum = atom->dvector[index_Atot_sum];
|
||||
ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
psi = atom->dvector[index_psi];
|
||||
psi_b = atom->dvector[index_psi_b];
|
||||
sigmaxx = atom->dvector[index_sigmaxx];
|
||||
sigmayy = atom->dvector[index_sigmayy];
|
||||
sigmazz = atom->dvector[index_sigmazz];
|
||||
history_setup_flag = atom->dvector[index_history_setup_flag];
|
||||
contacts = atom->dvector[index_contacts];
|
||||
adhesive_length = atom->dvector[index_adhesive_length];
|
||||
|
||||
end_of_step();
|
||||
}
|
||||
|
||||
void FixMDRradiusUpdate::pre_force(int)
|
||||
{
|
||||
|
||||
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
|
||||
class GranularModel* model;
|
||||
class GranularModel** models_list = pair->models_list;
|
||||
class GranSubModNormalMDR* norm_model = nullptr;
|
||||
for (int i = 0; i < pair->nmodels; i++) {
|
||||
model = models_list[i];
|
||||
if (model->normal_model->name == "mdr") norm_model = dynamic_cast<GranSubModNormalMDR *>(model->normal_model);
|
||||
}
|
||||
if (norm_model == nullptr) error->all(FLERR, "Did not find mdr model");
|
||||
//model = models_list[0];
|
||||
//class GranSubModNormalMDR* norm_model = dynamic_cast<GranSubModNormalMDR *>(model->normal_model);
|
||||
|
||||
//std::cout << "Preforce was called radius update" << std::endl;
|
||||
|
||||
// assign correct value to initially non-zero MDR particle history variables
|
||||
//int tmp1, tmp2;
|
||||
//int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
|
||||
//int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
|
||||
//int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
|
||||
//int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
|
||||
//int index_psi = atom->find_custom("psi",tmp1,tmp2);
|
||||
//int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
|
||||
//int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
|
||||
//int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
|
||||
//int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
|
||||
//int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
|
||||
//int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
|
||||
//int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
|
||||
//double * Ro = atom->dvector[index_Ro];
|
||||
//double * Vgeo = atom->dvector[index_Vgeo];
|
||||
//double * Velas = atom->dvector[index_Velas];
|
||||
//double * Atot = atom->dvector[index_Atot];
|
||||
//double * psi = atom->dvector[index_psi];
|
||||
//double * psi_b = atom->dvector[index_psi_b];
|
||||
//double * sigmaxx = atom->dvector[index_sigmaxx];
|
||||
//double * sigmayy = atom->dvector[index_sigmayy];
|
||||
//double * sigmazz = atom->dvector[index_sigmazz];
|
||||
//double * history_setup_flag = atom->dvector[index_history_setup_flag];
|
||||
//double * contacts = atom->dvector[index_contacts];
|
||||
//double * adhesive_length = atom->dvector[index_adhesive_length];
|
||||
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
//std::cout << "ntotal " << ntotal << ", nlocal " << atom->nlocal << ", nghost " << atom->nghost << std::endl;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (history_setup_flag[i] < 1e-16) {
|
||||
Ro[i] = radius[i];
|
||||
Vgeo[i] = 4.0/3.0*M_PI*pow(Ro[i],3.0);
|
||||
Velas[i] = 4.0/3.0*M_PI*pow(Ro[i],3.0);
|
||||
Atot[i] = 4.0*M_PI*pow(Ro[i],2.0);
|
||||
psi[i] = 1.0;
|
||||
psi_b[i] = norm_model->psi_b;
|
||||
history_setup_flag[i] = 1.0;
|
||||
}
|
||||
sigmaxx[i] = 0.0;
|
||||
sigmayy[i] = 0.0;
|
||||
sigmazz[i] = 0.0;
|
||||
contacts[i] = 0.0;
|
||||
adhesive_length[i] = 0.0;
|
||||
}
|
||||
|
||||
comm->forward_comm(this);
|
||||
|
||||
}
|
||||
|
||||
int FixMDRradiusUpdate::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
|
||||
{
|
||||
int m = 0;
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
buf[m++] = Ro[j]; // 1
|
||||
buf[m++] = Vgeo[j]; // 2
|
||||
buf[m++] = Velas[j]; // 3
|
||||
buf[m++] = Vcaps[j]; // 4
|
||||
buf[m++] = eps_bar[j]; // 5
|
||||
buf[m++] = dRnumerator[j]; // 6
|
||||
buf[m++] = dRdenominator[j]; // 7
|
||||
buf[m++] = Acon0[j]; // 8
|
||||
buf[m++] = Acon1[j]; // 9
|
||||
buf[m++] = Atot[j]; // 10
|
||||
buf[m++] = Atot_sum[j]; // 11
|
||||
buf[m++] = ddelta_bar[j]; // 12
|
||||
buf[m++] = psi[j]; // 13
|
||||
buf[m++] = psi_b[j]; // 14
|
||||
buf[m++] = sigmaxx[j]; // 15
|
||||
buf[m++] = sigmayy[j]; // 16
|
||||
buf[m++] = sigmazz[j]; // 17
|
||||
buf[m++] = history_setup_flag[j]; // 18
|
||||
buf[m++] = contacts[j]; // 19
|
||||
buf[m++] = adhesive_length[j]; // 20
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int m = 0;
|
||||
int last = first + n;
|
||||
for (int i = first; i < last; i++) {
|
||||
Ro[i] = buf[m++]; // 1
|
||||
Vgeo[i] = buf[m++]; // 2
|
||||
Velas[i] = buf[m++]; // 3
|
||||
Vcaps[i] = buf[m++]; // 4
|
||||
eps_bar[i] = buf[m++]; // 5
|
||||
dRnumerator[i] = buf[m++]; // 6
|
||||
dRdenominator[i] = buf[m++]; // 7
|
||||
Acon0[i] = buf[m++]; // 8
|
||||
Acon1[i] = buf[m++]; // 9
|
||||
Atot[i] = buf[m++]; // 10
|
||||
Atot_sum[i] = buf[m++]; // 11
|
||||
ddelta_bar[i] = buf[m++]; // 12
|
||||
psi[i] = buf[m++]; // 13
|
||||
psi_b[i] = buf[m++]; // 14
|
||||
sigmaxx[i] = buf[m++]; // 15
|
||||
sigmayy[i] = buf[m++]; // 16
|
||||
sigmazz[i] = buf[m++]; // 17
|
||||
history_setup_flag[i] = buf[m++]; // 18
|
||||
contacts[i] = buf[m++]; // 19
|
||||
adhesive_length[i] = buf[m++]; // 20
|
||||
}
|
||||
}
|
||||
|
||||
void FixMDRradiusUpdate::end_of_step()
|
||||
{
|
||||
|
||||
//comm->forward_comm(this);
|
||||
|
||||
//std::cout << "end_of_step() was called radius update" << std::endl;
|
||||
|
||||
// update the apparent radius of every particle
|
||||
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
double sigmaxx_sum = 0.0;
|
||||
double sigmayy_sum = 0.0;
|
||||
double sigmazz_sum = 0.0;
|
||||
double Vparticles = 0.0;
|
||||
|
||||
//std::cout << "New step" << std::endl;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
|
||||
const double R = radius[i];
|
||||
Atot[i] = 4.0*M_PI*pow(R,2.0) + Atot_sum[i];
|
||||
|
||||
const double Vo = 4.0/3.0*M_PI*pow(Ro[i],3.0);
|
||||
const double Vgeoi = 4.0/3.0*M_PI*pow(R,3.0) - Vcaps[i];
|
||||
Vgeo[i] = std::min(Vgeoi,Vo);
|
||||
//(Vgeoi < Vo) ? Vgeo[i] = Vgeoi : Vgeo[i] = Vo;
|
||||
|
||||
const double Afree = Atot[i] - Acon1[i];
|
||||
psi[i] = Afree/Atot[i];
|
||||
|
||||
//if (psi[i] < 0.08) {
|
||||
// std::cout << "psi is: " << psi[i] << std::endl;
|
||||
//}
|
||||
|
||||
const double dR = std::max(dRnumerator[i]/(dRdenominator[i] - 4.0*M_PI*pow(R,2.0)),0.0);
|
||||
if (psi_b[i] < psi[i]) {
|
||||
if ((radius[i] + dR) < (1.5*Ro[i])) radius[i] += dR;
|
||||
//radius[i] += dR;
|
||||
}
|
||||
|
||||
//if (dR > 1e-7) {
|
||||
// std::cout << "big dR change" << std::endl;
|
||||
//}
|
||||
//if (atom->tag[i] == 9){
|
||||
// std::cout << i << ", radius: " << radius[i] << ", dR: " << dR << ", dRnum: " << dRnumerator[i] << ", dRdenom " << dRdenominator[i] << ", dRdem_full " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl;
|
||||
//}
|
||||
|
||||
|
||||
//std::cout << i << ", " << radius[i] << " | " << psi_b[i] << ", " << psi[i] << " | " << Acon1[i] << ", " << Atot[i] << std::endl;
|
||||
|
||||
|
||||
Velas[i] = Vo*(1.0 + eps_bar[i]);
|
||||
Vcaps[i] = 0.0;
|
||||
eps_bar[i] = 0.0;
|
||||
dRnumerator[i] = 0.0;
|
||||
dRdenominator[i] = 0.0;
|
||||
Acon0[i] = Acon1[i];
|
||||
Acon1[i] = 0.0;
|
||||
//std::cout << "Acon reset: " << Acon0[i] << ", " << Acon1[i] << std::endl;
|
||||
Atot_sum[i] = 0.0;
|
||||
ddelta_bar[i] = 0.0;
|
||||
//adhesive_length[i] = adhesive_length[i]/contacts[i]; // convert adhesive length to average aAdh for each particle
|
||||
|
||||
|
||||
//if (lmp->update->ntimestep == 487500) {
|
||||
// double **x = atom->x;
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/atom_stresses_at_mid_compression.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
|
||||
// rowDataStream << sigmaxx[i] << ", " << sigmayy[i] << ", " << sigmazz[i] << ", " << Velas[i] << ", " << x[i][2];
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
//
|
||||
//if (lmp->update->ntimestep == 595750) {
|
||||
// double **x = atom->x;
|
||||
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/atom_stresses_at_max_compression.csv");
|
||||
// std::stringstream rowDataStream;
|
||||
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
|
||||
// rowDataStream << sigmaxx[i] << ", " << sigmayy[i] << ", " << sigmazz[i] << ", " << Velas[i] << ", " << x[i][2];
|
||||
// std::string rowData = rowDataStream.str();
|
||||
// csvWriter.writeRow(rowData);
|
||||
//}
|
||||
|
||||
|
||||
//sigmaxx_sum += sigmaxx[i];
|
||||
//sigmayy_sum += sigmayy[i];
|
||||
//sigmazz_sum += sigmazz[i];
|
||||
//Vparticles += Velas[i];
|
||||
}
|
||||
|
||||
//double sigmaxx_avg = sigmaxx_sum/nlocal;
|
||||
//double sigmayy_avg = sigmayy_sum/nlocal;
|
||||
//double sigmazz_avg = sigmazz_sum/nlocal;
|
||||
|
||||
comm->forward_comm(this);
|
||||
}
|
||||
|
||||
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_PI*pow(R,2.0) << std::endl;
|
||||
//std::cout << "Fix radius update setup has been entered !!!" << std::endl;
|
||||
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;
|
||||
66
src/fix_mdr_radius_update.h
Normal file
66
src/fix_mdr_radius_update.h
Normal file
@ -0,0 +1,66 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(mdr/radius/update,FixMDRradiusUpdate);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_MDR_RADIUS_UPDATE_H
|
||||
#define LMP_FIX_MDR_RADIUS_UPDATE_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixMDRradiusUpdate : public Fix {
|
||||
public:
|
||||
double * Ro;
|
||||
double * Vgeo;
|
||||
double * Velas;
|
||||
double * Vcaps;
|
||||
double * eps_bar;
|
||||
double * dRnumerator;
|
||||
double * dRdenominator;
|
||||
double * Acon0;
|
||||
double * Acon1;
|
||||
double * Atot;
|
||||
double * Atot_sum;
|
||||
double * ddelta_bar;
|
||||
double * psi;
|
||||
double * psi_b;
|
||||
double * sigmaxx;
|
||||
double * sigmayy;
|
||||
double * sigmazz;
|
||||
double * history_setup_flag;
|
||||
double * contacts;
|
||||
double * adhesive_length;
|
||||
|
||||
FixMDRradiusUpdate(class LAMMPS *, int, char **);
|
||||
int setmask() override;
|
||||
void setup(int) override;
|
||||
void setup_pre_force(int) override;
|
||||
void pre_force(int) override;
|
||||
void end_of_step() override; // FOR MDR
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
@ -77,7 +77,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
|
||||
SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
|
||||
RANDOM,NORMAL,CEIL,FLOOR,ROUND,TERNARY,
|
||||
RAMP,STAGGER,LOGFREQ,LOGFREQ2,LOGFREQ3,STRIDE,STRIDE2,
|
||||
VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,
|
||||
VDISPLACE,SWIGGLE,CWIGGLE,SIGN,GMASK,RMASK,
|
||||
GRMASK,IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,IS_FILE,EXTRACT_SETTING,
|
||||
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY};
|
||||
|
||||
@ -2575,7 +2575,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
|
||||
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),ternary(x,y,z),
|
||||
ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
|
||||
logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),swiggle(x,y,z),
|
||||
cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
|
||||
cwiggle(x,y,z),sign(x),gmask(x),rmask(x),grmask(x,y)
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
double Variable::collapse_tree(Tree *tree)
|
||||
@ -3125,6 +3125,14 @@ double Variable::collapse_tree(Tree *tree)
|
||||
return tree->value;
|
||||
}
|
||||
|
||||
if (tree->type == SIGN) {
|
||||
arg1 = collapse_tree(tree->first);
|
||||
if (tree->first->type != VALUE) return 0.0;
|
||||
tree->type = VALUE;
|
||||
tree->value = (arg1 >= 0.0) ? 1.0 : -1.0; // sign(arg1);
|
||||
return tree->value;
|
||||
}
|
||||
|
||||
// mask functions do not become a single collapsed value
|
||||
|
||||
if (tree->type == GMASK) return 0.0;
|
||||
@ -3143,7 +3151,7 @@ double Variable::collapse_tree(Tree *tree)
|
||||
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),ternary(x,y,z),
|
||||
ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
|
||||
logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),
|
||||
swiggle(x,y,z),cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
|
||||
swiggle(x,y,z),cwiggle(x,y,z),sign(x),gmask(x),rmask(x),grmask(x,y)
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
double Variable::eval_tree(Tree *tree, int i)
|
||||
@ -3459,6 +3467,9 @@ double Variable::eval_tree(Tree *tree, int i)
|
||||
return arg;
|
||||
}
|
||||
|
||||
if (tree->type == SIGN)
|
||||
return (eval_tree(tree->first,i) >= 0.0) ? 1.0 : -1.0; // sign(eval_tree(tree->first,i));
|
||||
|
||||
if (tree->type == GMASK) {
|
||||
if (atom->mask[i] & tree->ivalue) return 1.0;
|
||||
else return 0.0;
|
||||
@ -3650,7 +3661,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree
|
||||
strcmp(word,"logfreq") != 0 && strcmp(word,"logfreq2") != 0 &&
|
||||
strcmp(word,"logfreq3") != 0 && strcmp(word,"stride") != 0 &&
|
||||
strcmp(word,"stride2") != 0 && strcmp(word,"vdisplace") != 0 &&
|
||||
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0)
|
||||
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0)
|
||||
return 0;
|
||||
|
||||
// parse contents for comma-separated args
|
||||
@ -4045,6 +4056,11 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree
|
||||
double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
|
||||
argstack[nargstack++] = value;
|
||||
}
|
||||
} else if (strcmp(word,"sign") == 0) {
|
||||
if (narg != 1)
|
||||
print_var_error(FLERR,"Invalid math function in variable formula",ivar);
|
||||
if (tree) newtree->type = SIGN;
|
||||
else argstack[nargstack++] = (value1 >= 0.0) ? 1.0 : -1.0; // sign(value1);
|
||||
}
|
||||
|
||||
// delete stored args
|
||||
|
||||
@ -282,7 +282,7 @@ void Verlet::run(int n)
|
||||
}
|
||||
timer->stamp();
|
||||
comm->exchange();
|
||||
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
|
||||
//if (sortflag && ntimestep >= atom->nextsort) atom->sort();
|
||||
comm->borders();
|
||||
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
|
||||
timer->stamp(Timer::COMM);
|
||||
|
||||
Reference in New Issue
Block a user