commit e2df704fb21f7b57fc8060bd864671ac806569e5
parent 7170f28747ab13c1bb5994097d53fce776f5a129
Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
Date: Mon, 22 Nov 2021 10:15:25 +0100
fixed compaction and init, added layering.
Carpet mechanical property manipulation still needs to be tested.
Diffstat:
3 files changed, 118 insertions(+), 25 deletions(-)
diff --git a/compact_basin.jl b/compact_basin.jl
@@ -8,19 +8,23 @@ 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
+id = "simulation200" # id of simulation to load
N = 20e3 # amount of stress to be applied
-t_comp = 0.2 #compaction max duration [s]
+t_comp = 5.0 # compaction max duration [s]
sim = Granular.readSimulation("$(id)/init.jld2")
carpet = Granular.readSimulation("$(id)/carpet.jld2")
SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
+#mkpath("$(id)/compaction-N$(N)Pa")
-sim.id = "compaction-N$(N)Pa_Grains$(SimSettings["ngrains"])"
+cd("$id")
+sim.id = "compaction-N$(N)Pa"
+#sim.id = "$(id)/compaction-N$(N)Pa"
SimSettings["N"] = N
+
+
Granular.zeroKinematics!(sim)
Granular.zeroKinematics!(carpet)
@@ -40,13 +44,13 @@ Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top,
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
@@ -61,10 +65,9 @@ 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
+ for i = 1:100 # run for a while before measuring the state of the top wall
Granular.run!(sim, single_step=true)
end
@@ -74,8 +77,6 @@ while sim.time < sim.time_total
end
-
-
defined_normal_stress = (ones(size(effective_normal_stress,1)))
*(Granular.getWallNormalStress(sim, stress_type="effective"))
@@ -93,4 +94,8 @@ PyPlot.ylabel("Normal stress [Pa]")
PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
PyPlot.clf()
+cd("..")
+
+JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
+
Granular.writeSimulation(sim,filename = "$(id)/comp.jld2")
diff --git a/init_basin.jl b/init_basin.jl
@@ -3,7 +3,7 @@ import JLD2
import PyPlot
import Dates
-t_start = Dates.now() # Save the start time, print the end time later.
+t_start = Dates.now() # Save the start time, print the end time later.
############# Initialization Settings #############
@@ -12,13 +12,13 @@ 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
+ngrains = 200 # 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
+nx = aspect_ratio*ny*(stacks+1) # number of grain columns
ny = Int(round(ny/(stacks+1)))
nx = Int(round(nx/(stacks+1)))
@@ -54,11 +54,11 @@ SimSettings["r_max"] = r_max
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
+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,
@@ -109,8 +109,8 @@ for i = 1:stacks
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
+ if y_top < grain.lin_pos[2]
+ global y_top = grain.lin_pos[2]
end
end
@@ -118,9 +118,9 @@ for i = 1:stacks
# 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
+ 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,
@@ -226,8 +226,7 @@ Granular.writeSimulation(carpet,
JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
-
-#print time elapsed
+# print time elapsed
t_now = Dates.now()
dur = Dates.canonicalize(t_now-t_start)
print("Time elapsed: ",dur)
diff --git a/layer_basin.jl b/layer_basin.jl
@@ -0,0 +1,89 @@
+import Granular
+import JLD2
+import PyPlot
+import Dates
+
+id = "simulation200" # id of simulation to load, just write the folder
+ # name here
+
+# Layer interface positions
+# defined as decimals
+# ex: [0.4,0.6,1] will give two boundaries at 0.4 from the bottom
+# and 0.6 from the bottom. ie 3 layers
+# include 1 in the end as roof.
+interfaces = [0,0.4,0.6,1]
+
+# mechanical properties for each layer
+youngs_modulus = [2e7,2e7,2e7] # elastic modulus
+poissons_ratio = [0.185,0.200,0.185] # shear stiffness ratio
+tensile_strength = [0.0,0.0,0.0] # strength of bonds between grains
+contact_dynamic_friction = [0.4,0.4,0.4] # friction between grains
+rotating = [true,true,true] # can grains rotate or not
+
+sim = Granular.readSimulation("$(id)/comp.jld2")
+carpet = Granular.readSimulation("$(id)/carpet.jld2")
+SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
+
+Granular.zeroKinematics!(sim) # end any movement
+Granular.zeroKinematics!(carpet) # end any movement
+
+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
+
+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
+
+h = y_top-y_bot
+
+for grain in sim.grains
+
+ for i = 2:size(interfaces,1)
+
+ if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1]
+
+ grain.youngs_modulus = youngs_modulus[i-1]
+ grain.poissons_ratio = poissons_ratio[i-1]
+ grain.tensile_strength = tensile_strength[i-1]
+ grain.contact_dynamic_friction = contact_dynamic_friction[i-1]
+ grain.rotating = rotating[i-1]
+
+ end
+
+ end
+
+end
+
+#set the mechanical settings for the carpet
+for grain in carpet.grains
+ grain.youngs_modulus = youngs_modulus[2]
+ grain.poissons_ratio = poissons_ratio[2]
+ grain.tensile_strength = tensile_strength[2]
+ grain.contact_dynamic_friction = contact_dynamic_friction[2]
+ grain.rotating = rotating[2]
+end
+
+"""
+cd("$id")
+sim.id = "layered"
+
+Granular.resetTime!(sim)
+Granular.setTotalTime!(sim,0.3)
+Granular.run!(sim)
+
+cd("..")
+"""
+
+Granular.writeSimulation(sim,
+ filename = "$(id)/layered.jld2")
+
+Granular.writeSimulation(carpet,
+ filename = "$(id)/carpet.jld2")
