# ========================================================================= # Twisting a pre-buckled naturally curved rod # ========================================================================= variable N equal 300 # Beam geometry (mm) variable R equal 1000/44.84 variable L equal 300 variable l0 equal $(v_L/v_N) variable d equal 3.1 # Cross-section properties variable A equal (PI/4)*v_d^2 variable I equal (PI/64)*v_d^4 variable J equal (PI/32)*v_d^4 # Elastic moduli (kPa) variable E equal 1300 variable G equal 433 # ========================================================================= # Initialization # ========================================================================= units lj dimension 3 boundary f f f processors * 1 1 atom_style hybrid bpm/sphere dipole special_bonds lj 0.0 0.0 0.0 coul 1.0 1.0 1.0 newton on off neighbor $(5*v_l0) bin comm_modify mode single cutoff $(5*v_l0) vel yes read_data curved-rod-N${N}.lam # ========================================================================= # Groups # ========================================================================= group l1 id 1 group l2 id 2 group r2 id $(v_N) group r1 id $(v_N+1) group left union l1 l2 group right union r1 r2 group ends union l1 r1 group free subtract all ends group mobile union free right set group left type 2 set group right type 3 # ========================================================================= # BPM bond style # ========================================================================= variable Kr equal v_E*v_A/v_l0 variable Ks equal 12*v_E*v_I/v_l0^3 variable Kt equal v_G*v_J/v_l0 variable Kb equal v_E*v_I/v_l0 bond_style bpm/rotational smooth no break no overlay/pair no bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} & -1 -1 -1 -1 & 1e-2 1e-2 1e-2 1e-2 # ========================================================================= # Granular contact (Hertz-Mindlin) # ========================================================================= # Made up values variable YoungsModulus equal 10*v_E variable PoissonsRatio equal 0.4 variable FrictionCoeff equal 0.1 variable CoeffRestitution equal 0.7 pair_style hybrid zero $(5*v_l0) granular limit_damping yes pair_coeff * * zero pair_coeff 2 3 granular hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) & tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) & damping coeff_restitution # ========================================================================= # Computes and output # ========================================================================= compute mytemp all temp/sphere compute kerot all erotate/sphere # Note: variable torque and compute torque/chunk only works for group of point particles compute ft l1 property/atom fx fy fz tqx tqy tqz compute f l1 reduce sum c_ft[1] c_ft[2] c_ft[3] compute t l1 reduce sum c_ft[4] c_ft[5] c_ft[6] variable fx equal c_f[1] variable tx equal c_t[1] compute xmax all reduce max x compute xmin all reduce min x variable dx equal "c_xmax-c_xmin" compute quat all property/atom quatw quati quatj quatk compute ymaxabs all reduce maxabs y compute zmaxabs all reduce maxabs z variable cycles equal 0 thermo 10000 thermo_style custom step ke c_kerot pe v_cycles v_dx c_ymaxabs c_zmaxabs v_fx v_tx thermo_modify temp mytemp norm no # ========================================================================= # Stage 0: Straighten by applying end moment M = kappa * EI # (no gravity, no contact) # ========================================================================= dump 1 all custom 100000 ./bend-N${N}.dump & id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*] variable timestep equal 1e-3 timestep ${timestep} variable nsteps equal 5e6 variable Mz equal v_E*v_I/v_R variable mz equal -v_Mz*ramp(0,1) fix 1 mobile nve/bpm/sphere update dipole fix 2 r1 addtorque/atom 0 0 v_mz fix 3 free gravity 0 vector 0 -1 0 fix 4 free langevin 0 0 1e3 12345 omega yes fix 5 all setforce NULL NULL 0 run $(v_nsteps) variable mz equal -v_Mz fix 4 free langevin 0 0 1e2 12345 omega yes run $(14*v_nsteps) unfix 1 unfix 2 unfix 3 unfix 4 unfix 5 undump 1 set group all dipole 0 1 0 # write_restart straightened-beam.restart # ========================================================================= # Stage 1: Axial compression (80 mm) # ========================================================================= balance 1.2 shift x 10 1.1 pair_coeff 1 1 granular hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) & tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) & damping coeff_restitution dump 1 all custom 100000 ./compress-N${N}.dump & id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*] variable timestep equal 1e-5 timestep ${timestep} variable nsteps equal 5e7 variable period equal v_timestep*v_nsteps variable vx equal -80/(v_timestep*v_nsteps) fix 1 free nve/bpm/sphere update dipole fix 2 r1 move linear ${vx} 0 0 fix 3 free gravity 9.81e-3 vector 0 -1 0 fix 4 free langevin 0 0 1e2 12345 omega yes run $(v_nsteps) # ========================================================================= # Stage 2: Relaxation # ========================================================================= unfix 2 fix 4 free langevin 0 0 1e3 12345 omega yes run $(4*v_nsteps) fix 4 free langevin 0 0 1e2 12345 omega yes run $(5*v_nsteps) # write_restart compressed-beam.restart undump 1 unfix 1 unfix 3 unfix 4 # ========================================================================= # Stage 3: Twist (9 full rotations, with wall) # ========================================================================= variable timestep equal 1e-4 timestep ${timestep} reset_timestep 0 balance 1.2 shift x 10 1.1 dump 1 all custom 1000000 ./twist-N${N}.dump & id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*] variable nturns equal 9 variable stepsperturn equal 1e8 variable period equal v_timestep*v_stepsperturn variable cycles equal step/v_stepsperturn variable theta equal 2*PI*v_cycles fix 1 free nve/bpm/sphere update dipole fix 2 r1 move rotate 0 0 0 -1 0 0 $(v_period) update dipole fix 3 free gravity 9.81e-3 vector 0 -1 0 fix 4 free langevin 0 0 1e2 12345 omega yes fix 5 free wall/gran granular & hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) & tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) & damping coeff_restitution & xplane 0 220 fix 6 all print 100000 & "${theta} $(c_ymaxabs) $(c_zmaxabs) ${tx}" & file twist-N${N}-theta-yzmaxabs-Tx.txt screen no & title "theta [rad], y_maxabs [mm], z_maxabs [mm], Tx [uJ]" variable i loop ${nturns} label loopstart run $(v_stepsperturn) next i jump SELF loopstart