shear.jl (10814B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import JLD2 5 import PyPlot 6 7 ################################################################################ 8 #### SIMULATION PARAMETERS # 9 ################################################################################ 10 let 11 12 # Common simulation identifier 13 id_prefix = "test0" 14 15 # Gravitational acceleration vector (cannot be zero; required for Step 1) 16 g = [0., -9.8] 17 18 # Grain package geometry during initialization 19 nx = 10 # Grains along x (horizontal) 20 ny = 50 # Grains along y (vertical) 21 22 # Grain-size parameters 23 r_min = 0.03 # Min. grain radius [m] 24 r_max = 0.1 # Max. grain radius [m] 25 gsd_type = "powerlaw" # "powerlaw" or "uniform" sizes between r_min and r_max 26 gsd_powerlaw_exponent = -1.8 # GSD power-law exponent 27 gsd_seed = 1 # Value to seed random-size generation 28 29 # Grain mechanical properties 30 youngs_modulus = 2e7 # Elastic modulus [Pa] 31 poissons_ratio = 0.185 # Shear-stiffness ratio [-] 32 tensile_strength = 0.0 # Inter-grain bond strength [Pa] 33 contact_dynamic_friction = 0.4 # Coulomb-frictional coefficient [-] 34 rotating = true # Allow grain rotation 35 36 # Normal stress for the consolidation and shear [Pa] 37 N = 20e3 38 39 # Shear velocity to apply to the top grains [m/s] 40 vel_shear = 0.5 41 42 # Simulation duration of individual steps [s] 43 t_init = 2.0 44 t_cons = 2.5 45 t_shear = 5.0 46 47 ################################################################################ 48 #### Step 1: Create a loose granular assemblage and let it settle at -y # 49 ################################################################################ 50 sim = Granular.createSimulation(id="$(id_prefix)-init") 51 52 Granular.regularPacking!(sim, [nx, ny], r_min, r_max, 53 size_distribution=gsd_type, 54 size_distribution_parameter=gsd_powerlaw_exponent, 55 seed=gsd_seed) 56 57 # Set grain mechanical properties 58 for grain in sim.grains 59 grain.youngs_modulus = youngs_modulus 60 grain.poissons_ratio = poissons_ratio 61 grain.tensile_strength = tensile_strength 62 grain.contact_dynamic_friction = contact_dynamic_friction 63 grain.rotating = rotating 64 end 65 66 # Create a grid for contact searching spanning the extent of the grains 67 Granular.fitGridToGrains!(sim, sim.ocean) 68 69 # Make the top and bottom boundaries impermeable, and the side boundaries 70 # periodic, which will come in handy during shear 71 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 72 verbose=false) 73 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west") 74 75 # Add gravitational acceleration to all grains and disable ocean-grid drag. 76 # Also add viscous energy dissipation between grains, which is disabled before 77 # consolidation and shear. 78 for grain in sim.grains 79 Granular.addBodyForce!(grain, grain.mass*g) 80 Granular.disableOceanDrag!(grain) 81 grain.contact_viscosity_normal = 1e4 # N/(m/s) 82 end 83 84 # Automatically set the computational time step based on grain sizes and 85 # properties 86 Granular.setTimeStep!(sim) 87 88 # Set the total simulation time for this step [s] 89 # This value may need tweaking if grain sizes or numbers are adjusted. 90 Granular.setTotalTime!(sim, t_init) 91 92 # Set the interval in model time between simulation files [s] 93 Granular.setOutputFileInterval!(sim, .02) 94 95 # Visualize the grain-size distribution 96 Granular.plotGrainSizeDistribution(sim) 97 98 # Start the simulation 99 Granular.run!(sim) 100 101 # Try to render the simulation if `pvpython` is installed on the system 102 Granular.render(sim, trim=false) 103 104 # Save the simulation state to disk in case we need to reuse the current state 105 Granular.writeSimulation(sim) 106 107 # Also copy the simulation in memory, in case we want to loop over different 108 # normal stresses below: 109 sim_init = deepcopy(sim) 110 111 112 ################################################################################ 113 #### Step 2: Consolidate the previous output under a constant normal stress # 114 ################################################################################ 115 116 # Rename the simulation so we don't overwrite output from the previous step 117 sim.id = "$(id_prefix)-cons-N$(N)Pa" 118 119 # Set all linear and rotational velocities to zero 120 Granular.zeroKinematics!(sim) 121 122 # Add a dynamic wall to the top which adds a normal stress downwards. The 123 # normal of this wall is downwards, and we place it at the top of the granular 124 # assemblage. Here, the inter-grain viscosity is also removed. 125 y_top = -Inf 126 for grain in sim.grains 127 grain.contact_viscosity_normal = 0. 128 if y_top < grain.lin_pos[2] + grain.contact_radius 129 y_top = grain.lin_pos[2] + grain.contact_radius 130 end 131 end 132 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top, 133 bc="normal stress", normal_stress=-N, 134 contact_viscosity_normal=1e3) 135 @info "Placing top wall at y=$y_top" 136 137 # Resize the grid to span the current state 138 Granular.fitGridToGrains!(sim, sim.ocean) 139 140 # Lock the grains at the very bottom so that the lower boundary is rough 141 y_bot = Inf 142 for grain in sim.grains 143 if y_bot > grain.lin_pos[2] - grain.contact_radius 144 y_bot = grain.lin_pos[2] - grain.contact_radius 145 end 146 end 147 fixed_thickness = 2. * r_max 148 for grain in sim.grains 149 if grain.lin_pos[2] <= fixed_thickness 150 grain.fixed = true # set x and y acceleration to zero 151 end 152 end 153 154 # Set current time to zero and reset output file counter 155 Granular.resetTime!(sim) 156 157 # Set the simulation time to run the consolidation for 158 Granular.setTotalTime!(sim, t_cons) 159 160 # Run the consolidation experiment, and monitor top wall position over time 161 time = Float64[] 162 compaction = Float64[] 163 effective_normal_stress = Float64[] 164 while sim.time < sim.time_total 165 166 for i=1:100 # run for 100 steps before measuring shear stress and dilation 167 Granular.run!(sim, single_step=true) 168 end 169 170 append!(time, sim.time) 171 append!(compaction, sim.walls[1].pos) 172 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 173 174 end 175 defined_normal_stress = ones(length(effective_normal_stress)) * 176 Granular.getWallNormalStress(sim, stress_type="effective") 177 PyPlot.subplot(211) 178 PyPlot.subplots_adjust(hspace=0.0) 179 ax1 = PyPlot.gca() 180 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 181 PyPlot.plot(time, compaction) 182 PyPlot.ylabel("Top wall height [m]") 183 PyPlot.subplot(212, sharex=ax1) 184 PyPlot.plot(time, defined_normal_stress) 185 PyPlot.plot(time, effective_normal_stress) 186 PyPlot.xlabel("Time [s]") 187 PyPlot.ylabel("Normal stress [Pa]") 188 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 189 PyPlot.clf() 190 191 # Try to render the simulation if `pvpython` is installed on the system 192 Granular.render(sim, trim=false) 193 194 # Save the simulation state to disk in case we need to reuse the consolidated 195 # state (e.g. different shear velocities below) 196 Granular.writeSimulation(sim) 197 198 # Also copy the simulation in memory, in case we want to loop over different 199 # normal stresses below: 200 sim_cons = deepcopy(sim) 201 202 203 ################################################################################ 204 #### Step 3: Shear the consolidated assemblage with a constant velocity # 205 ################################################################################ 206 207 # Rename the simulation so we don't overwrite output from the previous step 208 sim.id = "$(id_prefix)-shear-N$(N)Pa-vel_shear$(vel_shear)m-s" 209 210 # Set all linear and rotational velocities to zero 211 Granular.zeroKinematics!(sim) 212 213 # Set current time to zero and reset output file counter 214 Granular.resetTime!(sim) 215 216 # Set the simulation time to run the shear experiment for 217 Granular.setTotalTime!(sim, t_shear) 218 219 # Run the shear experiment 220 time = Float64[] 221 shear_stress = Float64[] 222 shear_strain = Float64[] 223 dilation = Float64[] 224 thickness_initial = sim.walls[1].pos - y_bot 225 x_min = +Inf 226 x_max = -Inf 227 for grain in sim.grains 228 if x_min > grain.lin_pos[1] - grain.contact_radius 229 x_min = grain.lin_pos[1] - grain.contact_radius 230 end 231 if x_max < grain.lin_pos[1] + grain.contact_radius 232 x_max = grain.lin_pos[1] + grain.contact_radius 233 end 234 end 235 surface_area = (x_max - x_min) 236 shear_force = 0. 237 while sim.time < sim.time_total 238 239 # Prescribe the shear velocity to the uppermost grains 240 for grain in sim.grains 241 if grain.lin_pos[2] >= sim.walls[1].pos - fixed_thickness 242 grain.fixed = true 243 grain.allow_y_acc = true 244 grain.lin_vel[1] = vel_shear 245 elseif grain.lin_pos[2] > fixed_thickness # do not free bottom grains 246 grain.fixed = false 247 end 248 end 249 250 for i=1:100 # run for 100 steps before measuring shear stress and dilation 251 Granular.run!(sim, single_step=true) 252 end 253 254 append!(time, sim.time) 255 256 # Determine the current shear stress 257 shear_force = 0. 258 for grain in sim.grains 259 if grain.fixed && grain.allow_y_acc 260 shear_force += -grain.force[1] 261 end 262 end 263 append!(shear_stress, shear_force/surface_area) 264 265 # Determine the current shear strain 266 append!(shear_strain, sim.time*vel_shear/thickness_initial) 267 268 # Determine the current dilation 269 append!(dilation, (sim.walls[1].pos - y_bot)/thickness_initial) 270 271 end 272 273 # Try to render the simulation if `pvpython` is installed on the system 274 Granular.render(sim, trim=false) 275 276 # Save the simulation state to disk in case we need to reuse the sheared state 277 Granular.writeSimulation(sim) 278 279 # Plot time vs. shear stress and dilation 280 PyPlot.subplot(211) 281 PyPlot.subplots_adjust(hspace=0.0) 282 ax1 = PyPlot.gca() 283 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 284 PyPlot.plot(time, shear_stress) 285 PyPlot.ylabel("Shear stress [Pa]") 286 PyPlot.subplot(212, sharex=ax1) 287 PyPlot.plot(time, dilation) 288 PyPlot.xlabel("Time [s]") 289 PyPlot.ylabel("Volumetric strain [-]") 290 PyPlot.savefig(sim.id * "-time_vs_shear-stress.pdf") 291 PyPlot.clf() 292 293 # Plot shear strain vs. shear stress and dilation 294 PyPlot.subplot(211) 295 PyPlot.subplots_adjust(hspace=0.0) 296 ax1 = PyPlot.gca() 297 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 298 PyPlot.plot(shear_strain, shear_stress) 299 PyPlot.ylabel("Shear stress [Pa]") 300 PyPlot.subplot(212, sharex=ax1) 301 PyPlot.plot(shear_strain, dilation) 302 PyPlot.xlabel("Shear strain [-]") 303 PyPlot.ylabel("Volumetric strain [-]") 304 PyPlot.savefig(sim.id * "-shear-strain_vs_shear-stress.pdf") 305 PyPlot.clf() 306 307 # Plot time vs. shear strain (boring when the shear experiment is rate controlled) 308 PyPlot.plot(time, shear_strain) 309 PyPlot.xlabel("Time [s]") 310 PyPlot.ylabel("Shear strain [-]") 311 PyPlot.savefig(sim.id * "-time_vs_shear-strain.pdf") 312 PyPlot.clf() 313 end