commit 3cef6c25316427fe9461d3f64ebded471c61de96
parent 125e84ccadfa1ace469d3fe5e2b64e9969770e28
Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
Date: Tue, 23 Nov 2021 17:25:24 +0100
minor changes.
Large additions to deform_basin.jl
Diffstat:
4 files changed, 117 insertions(+), 30 deletions(-)
diff --git a/compact_basin.jl b/compact_basin.jl
@@ -10,7 +10,7 @@ t_start = Dates.now() # Save the start time, print the end time later.
id = "simulation1000" # id of simulation to load
N = 20e3 # amount of stress to be applied
-t_comp = 5.0 # compaction max duration [s]
+t_comp = 3.0 # compaction max duration [s]
sim = Granular.readSimulation("$(id)/init.jld2")
carpet = Granular.readSimulation("$(id)/carpet.jld2")
@@ -24,7 +24,6 @@ sim.id = "compaction-N$(N)Pa"
SimSettings["N"] = N
-
Granular.zeroKinematics!(sim)
Granular.zeroKinematics!(carpet)
diff --git a/deform_basin.jl b/deform_basin.jl
@@ -10,12 +10,68 @@ id = "simulation1000" # folder name of simulation
hw_ratio = 0.2 # height/width ratio of indenter
grain_radius = 0.05 # grain radius of grains in indenter
-deformation_type = # "diapir" or "inversion"
+deformation_type = "diapir" # "diapir" or "inversion"
+ # diapir will only introduce an indenter while
+ # inversion will also add moving east/west walls
+ # that follow the contraction of the carpet
t_start = Dates.now()
-sim = Granular.readSimulation("$(id)/comp.jld2")
-carpet = Granular.readSimulation("$(id)/carpet.jld2")
+sim = Granular.readSimulation("$(id)/layered.jld2")
SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
-temp_indent = createSimulation("id=temp_indent")
+# Add Indenter
+temp_indent = Granular.createSimulation("id=temp_indent")
+
+left_edge = round(sim.ocean.origo[1],digits=2)
+length = round(sim.ocean.L[1],digits=2)
+
+width = length/3
+hw_ratio = 0.2
+init_vertex_pos = [(length+left_edge)/2,-0.2]
+grain_radius = 0.05
+
+vertex_x = init_vertex_pos[1]
+vertex_y = width*hw_ratio*sin((pi/width)*vertex_x)
+
+
+for i = 0:grain_radius*2:width#manipulate the ocean grid
+
+ x_pos = i
+
+ y_pos = width*hw_ratio*sin(pi/width*x_pos)
+
+ Granular.addGrainCylindrical!(temp_indent,
+ [x_pos+init_vertex_pos[1]-width/2,y_pos+vertex_y+init_vertex_pos[2]],
+ grain_radius,
+ 0.1,
+ fixed = true,
+ lin_vel = [0.0,0.5])
+end
+
+append!(sim.grains,temp_indent.grains)
+
+Granular.fitGridToGrains!(sim,
+ sim.ocean,
+ north_padding = 3.0,
+ verbose=false)
+
+sim.time_iteration = 0
+sim.time = 0.0
+sim.file_time_since_output_file = 0.
+Granular.setTotalTime!(sim,2.0)
+Granular.setTimeStep!(sim)
+Granular.setOutputFileInterval!(sim, .01)
+
+
+cd("$id")
+sim.id = "deformed"
+
+Granular.resetTime!(sim)
+Granular.setTotalTime!(sim,2.0)
+Granular.run!(sim)
+
+cd("..")
+
+Granular.writeSimulation(sim,
+ filename = "$(id)/deformed.jld2")
diff --git a/init_basin.jl b/init_basin.jl
@@ -12,7 +12,7 @@ t_stack = 0.8 # 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
+ngrains = 1000 # total number of grains
aspect_ratio = 4 # should be x times as wide as it is tall
mkpath("simulation$(ngrains)")
@@ -186,6 +186,7 @@ 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
+# color = 1 is used as a flag for the grains in the carpet
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
@@ -199,7 +200,8 @@ for i = left_edge+(bot_r/2):bot_r*1.99:left_edge+length
shear_strength = Inf,
contact_stiffness_normal = Inf,
contact_stiffness_tangential = Inf,
- fixed = true)
+ fixed = true,
+ color = 1)
end
Granular.findContactsAllToAll!(carpet) # find the grain contacts
diff --git a/layer_basin.jl b/layer_basin.jl
@@ -15,17 +15,17 @@ 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
+poissons_ratio = [0.185,0.185,0.185] # shear stiffness ratio
+tensile_strength = [0.3,0.05,0.3] # 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
+color = [0,0,0]
-#mechanical properties for carpet
-carpet_youngs_modulus = 2e7 # elastic modulus
-carpet_poissons_ratio = 0.185 # shear stiffness ratio
-carpet_tensile_strength = Inf # strength of bonds between grains
-carpet_contact_dynamic_friction = 0.4 # friction between grains
-carpet_rotating = true # can grains rotate or not
+carpet_youngs_modulus = 2e7
+carpet_poissons_ratio = 0.185
+carpet_tensile_strength = Inf
+carpet_contact_dynamic_friction = 0.4
+carpet_rotating = true
sim = Granular.readSimulation("$(id)/comp.jld2")
carpet = Granular.readSimulation("$(id)/carpet.jld2")
@@ -49,46 +49,76 @@ for grain in sim.grains
end
end
-h = y_top-y_bot
+# quick fix to make the color = 1 flag for grains belonging to the carpet.
+# this should be done in the newer versions of init_basin.jl instead
+"""
+for grain in sim.grains
+ if grain.lin_pos[2] == -0.05
+ grain.color = 1
+ end
+end
+"""
+
+h = y_top-y_bot #depth of basin
+
+interfaces *= h
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]
+ if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1] && grain.color != 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]
-
+ grain.color = color[i-1]
+ elseif grain.color == 1
+ grain.youngs_modulus = carpet_youngs_modulus
+ grain.poissons_ratio = carpet_poissons_ratio
+ grain.tensile_strength = carpet_tensile_strength
+ grain.contact_dynamic_friction = carpet_contact_dynamic_friction
+ grain.rotating = carpet_rotating
end
-
end
-
end
-#set the mechanical settings for the carpet
-for grain in carpet.grains
- grain.youngs_modulus = carpet_youngs_modulus
- grain.poissons_ratio = carpet_poissons_ratio
- grain.tensile_strength = carpet_tensile_strength
- grain.contact_dynamic_friction = carpet_contact_dynamic_friction
- grain.rotating = carpet_rotating
+# Create the contacs between grains by expanding all grains by a small amount
+# then search and establish contacts and then reduce the size of the grains again
+
+size_increasing_factor = 1.10 # factor by which contact radius should be increased
+ # to search for contacts
+increase_array = []
+
+#increase the contact radius
+for grain in sim.grains
+ if grain.color == 0
+ contact_radius_increase = (grain.contact_radius*size_increasing_factor)-grain.contact_radius
+ grain.contact_radius += contact_radius_increase
+ append!(increase_array,contact_radius_increase)
+ elseif grain.color == 1
+ append!(increase_array,0)
+ end
end
+Granular.findContactsAllToAll!(carpet) # find the grain contacts
+
+#reduce the contact radius again
+for i = 1:size(sim.grains,1)
+ sim.grains[i].contact_radius -= increase_array[i]
+end
cd("$id")
sim.id = "layered"
Granular.resetTime!(sim)
-Granular.setTotalTime!(sim,0.05)
+Granular.setTotalTime!(sim,0.2)
Granular.run!(sim)
cd("..")
-
Granular.writeSimulation(sim,
filename = "$(id)/layered.jld2")
