tinitial commit of init and compact files - granular-basin - tectonic deformation experiments with Granular.jl
(HTM) git clone git://src.adamsgaard.dk/granular-basin
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
---
(DIR) commit 7170f28747ab13c1bb5994097d53fce776f5a129
(DIR) parent 8a5aed9e374b8484f1695b47dbbeacb1496b0a59
(HTM) Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
Date: Wed, 17 Nov 2021 16:05:01 +0100
initial commit of init and compact files
Diffstat:
A compact_basin.jl | 96 +++++++++++++++++++++++++++++++
A init_basin.jl | 233 +++++++++++++++++++++++++++++++
2 files changed, 329 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/compact_basin.jl b/compact_basin.jl
t@@ -0,0 +1,96 @@
+import Granular
+import JLD2
+import PyPlot
+import Dates
+
+t_start = Dates.now() # Save the start time, print the end time later.
+
+# lav en lille test? se om dit appendede carpet stadig er forbundet til hoved-
+# simulationsobjektet
+
+
+id = "simulation500" # id of simulation to load
+N = 20e3 # amount of stress to be applied
+t_comp = 0.2 #compaction max duration [s]
+
+sim = Granular.readSimulation("$(id)/init.jld2")
+carpet = Granular.readSimulation("$(id)/carpet.jld2")
+SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
+
+
+sim.id = "compaction-N$(N)Pa_Grains$(SimSettings["ngrains"])"
+SimSettings["N"] = N
+
+Granular.zeroKinematics!(sim)
+
+Granular.zeroKinematics!(carpet)
+
+y_top = -Inf
+for grain in sim.grains
+ grain.contact_viscosity_normal = 0
+ if y_top < grain.lin_pos[2] + grain.contact_radius
+ global y_top = grain.lin_pos[2] + grain.contact_radius
+ end
+end
+
+Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top,
+ bc="normal stress",
+ normal_stress=-N,
+ contact_viscosity_normal=1e3)
+
+Granular.fitGridToGrains!(sim,sim.ocean)
+
+
+y_bot = Inf
+for grain in sim.grains
+ if y_bot > grain.lin_pos[2] - grain.contact_radius
+ global y_bot = grain.lin_pos[2] - grain.contact_radius
+ end
+end
+fixed_thickness = 2. * SimSettings["r_max"]
+for grain in sim.grains
+ if grain.lin_pos[2] <= fixed_thickness
+ grain.fixed = true # set x and y acceleration to zero
+ end
+end
+
+Granular.resetTime!(sim)
+Granular.setTotalTime!(sim,t_comp)
+
+time = Float64[]
+compaction = Float64[]
+effective_normal_stress = Float64[]
+
+
+while sim.time < sim.time_total
+
+ for i = 1:100 #run for a while before measuring the state of the top wall
+ Granular.run!(sim, single_step=true)
+ end
+
+ append!(time, sim.time)
+ append!(compaction, sim.walls[1].pos)
+ append!(effective_normal_stress, Granular.getWallNormalStress(sim))
+
+end
+
+
+
+defined_normal_stress = (ones(size(effective_normal_stress,1)))
+ *(Granular.getWallNormalStress(sim, stress_type="effective"))
+
+PyPlot.subplot(211)
+PyPlot.subplots_adjust(hspace=0.0)
+ax1 = PyPlot.gca()
+PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels
+PyPlot.plot(time, compaction)
+PyPlot.ylabel("Top wall height [m]")
+PyPlot.subplot(212, sharex=ax1)
+PyPlot.plot(time, defined_normal_stress)
+PyPlot.plot(time, effective_normal_stress)
+PyPlot.xlabel("Time [s]")
+PyPlot.ylabel("Normal stress [Pa]")
+PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
+PyPlot.clf()
+
+Granular.writeSimulation(sim,filename = "$(id)/comp.jld2")
(DIR) diff --git a/init_basin.jl b/init_basin.jl
t@@ -0,0 +1,233 @@
+import Granular
+import JLD2
+import PyPlot
+import Dates
+
+t_start = Dates.now() # Save the start time, print the end time later.
+
+############# Initialization Settings #############
+
+t_init = 1.0 # duration of initialization [s]
+t_stack = 0.5 # duration for each stack to settle [s]
+
+g = [0.,-9.8] # vector for direction and magnitude of gravitational acceleration of grains
+
+ngrains = 40000 # total number of grains
+aspect_ratio = 2 #should be x times as wide as it is tall
+
+stacks = 2 # number of duplicate stacks on top of the initial grains
+
+ny = sqrt((ngrains)/aspect_ratio) # number of grain rows
+nx = aspect_ratio*ny*(stacks+1) # number of grain columns
+
+ny = Int(round(ny/(stacks+1)))
+nx = Int(round(nx/(stacks+1)))
+
+
+r_min = 0.05 # minimum radius of grains
+r_max = r_min*sqrt(2) # max radius of grains, double the area of minimum
+
+# Grain-size distribution parameters
+gsd_type = "powerlaw" # type grain-size distribution
+gsd_powerlaw_exponent = -1.8 # exponent if powerlaw is used for grain-size distribution
+gsd_seed = 3 # random seed for the distribution of grain-sizes
+
+# Initial mechanical properties of grains
+youngs_modulus = 2e7 # elastic modulus
+poissons_ratio = 0.185 # shear stiffness ratio
+tensile_strength = 0.0 # strength of bonds between grains
+contact_dynamic_friction = 0.4 # friction between grains
+rotating = true # can grains rotate or not
+
+# Save some settings in a dictionary
+SimSettings = Dict()
+SimSettings["nx"] = nx
+SimSettings["ny"] = ny
+SimSettings["stacks"] = stacks
+SimSettings["ngrains"] = ngrains
+SimSettings["r_min"] = r_min
+SimSettings["r_max"] = r_max
+
+
+
+############# Initialize simulation and grains #############
+
+sim = Granular.createSimulation(id="simulation$(ngrains)")
+
+Granular.regularPacking!(sim, #simulation object
+ [nx, ny], #number of grains along x and y axis
+ r_min, #smallest grain size
+ r_max, #largest grain size
+ tiling="triangular", #how the grains are tiled
+ origo = [0.0,0.0],
+ size_distribution=gsd_type,
+ size_distribution_parameter=gsd_powerlaw_exponent,
+ seed=gsd_seed)
+
+
+# set the indicated mechanical parameters for all grains
+for grain in sim.grains
+ grain.youngs_modulus = youngs_modulus
+ grain.poissons_ratio = poissons_ratio
+ grain.tensile_strength = tensile_strength
+ grain.contact_dynamic_friction = contact_dynamic_friction
+ grain.rotating = rotating
+end
+
+# fit the ocean grid to the grains
+Granular.fitGridToGrains!(sim, sim.ocean)
+
+# set the boundary conditions
+Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
+ verbose=false)
+Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
+ verbose=false)
+
+# add gravity to the grains and remove ocean drag
+for grain in sim.grains
+ Granular.addBodyForce!(grain, grain.mass*g)
+ Granular.disableOceanDrag!(grain)
+ grain.contact_viscosity_normal = 1e4 # N/(m/s)
+end
+
+# run the simulation
+Granular.setTimeStep!(sim) # set appropriate time steps
+Granular.setTotalTime!(sim, t_init) # set total time
+Granular.setOutputFileInterval!(sim, .01) # how often vtu files should be outputted
+
+Granular.run!(sim)
+
+
+
+############# Stack the initialized grains #############
+
+temp = deepcopy(sim)
+
+for i = 1:stacks
+
+ # find the position of the uppermost grain
+ global y_top = -Inf
+
+ for grain in sim.grains
+ if y_top < grain.lin_pos[2] #+ grain.contact_radius
+ global y_top = grain.lin_pos[2] #+ grain.contact_radius
+ end
+ end
+
+
+ # add duplicate grains above the initialized grains
+ for grain in temp.grains
+
+ lin_pos_lifted = [0.0,0.0] # preallocation of position of a 'lifted' grain
+ lin_pos_lifted[1] = grain.lin_pos[1] # x position of duplicate grain
+ lin_pos_lifted[2] = grain.lin_pos[2] + y_top + r_max # y-position of duplicate grain
+
+
+ Granular.addGrainCylindrical!(sim,
+ lin_pos_lifted,
+ grain.contact_radius,
+ grain.thickness,
+ youngs_modulus = grain.youngs_modulus,
+ poissons_ratio = grain.poissons_ratio,
+ tensile_strength = grain.tensile_strength,
+ contact_dynamic_friction = grain.contact_dynamic_friction,
+ verbose = false)
+
+ end
+
+ Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
+ verbose=false)
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west",
+ verbose=false)
+ Granular.setTotalTime!(sim,t_stack)
+ Granular.setTimeStep!(sim)
+ Granular.setOutputFileInterval!(sim, .01)
+
+ for grain in sim.grains
+ Granular.addBodyForce!(grain, grain.mass*g)
+ Granular.disableOceanDrag!(grain)
+ end
+
+
+ # instead of Granular.resetTime!(sim), the time parameters are reset manually, in order
+ # to allow for the continuation of vtu-file index from earlier
+ sim.time_iteration = 0
+ sim.time = 0.0
+ sim.file_time_since_output_file = 0.
+
+ Granular.setTotalTime!(sim, t_stack)
+ Granular.setTimeStep!(sim)
+ Granular.setOutputFileInterval!(sim, .01)
+
+ Granular.run!(sim) # let the duplicated grains settle.
+
+end
+
+
+
+############# Lay a carpet #############
+
+carpet = Granular.createSimulation(id="init_carpet") # new simulation object for the carpet
+
+bot_r = r_min # radius of carpet grains
+
+left_edge = round(sim.ocean.origo[1],digits=2) # west edge of the carpet
+length = round(sim.ocean.L[1],digits=2) # width of the carpet
+right_edge = left_edge+length # east edge of the carpet
+
+
+# Now loop over the carpet grain positions, the loop will create grains that overlap slightly
+# in order to create the bonds needed
+for i = left_edge+(bot_r/2):bot_r*1.99:left_edge+length
+
+ bot_pos = [i,round(sim.ocean.origo[2]-bot_r,digits=2)] # position of grain
+
+ Granular.addGrainCylindrical!(carpet,
+ bot_pos,
+ bot_r,
+ 0.1,
+ verbose = false,
+ tensile_strength = Inf,
+ shear_strength = Inf,
+ contact_stiffness_normal = Inf,
+ contact_stiffness_tangential = Inf,
+ fixed = true)
+end
+
+Granular.findContactsAllToAll!(carpet) # find the grain contacts
+
+append!(sim.grains,carpet.grains) # add the carpet grains to the main simulation object
+# since the assignment will point to the carpet object, changes made to the carpet
+# object will appear in the main simulation object
+
+
+Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) # fit the ocean to the added grains
+
+# run the simulation shortly, to let the stacked grains settle on the carpet
+sim.time_iteration = 0
+sim.time = 0.0
+sim.file_time_since_output_file = 0.
+Granular.setTotalTime!(sim, 0.5)
+Granular.setTimeStep!(sim)
+Granular.setOutputFileInterval!(sim, .01)
+
+Granular.run!(sim)
+
+
+# save the simulation and the carpet objects
+
+Granular.writeSimulation(sim,
+ filename = "simulation$(ngrains)/init.jld2")
+
+
+Granular.writeSimulation(carpet,
+ filename = "simulation$(ngrains)/carpet.jld2")
+
+JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
+
+
+#print time elapsed
+t_now = Dates.now()
+dur = Dates.canonicalize(t_now-t_start)
+print("Time elapsed: ",dur)