tlayer_basin.jl - 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
---
tlayer_basin.jl (3778B)
---
1 import Granular
2 import JLD2
3 import PyPlot
4 import Dates
5
6 id = "simulation40000" # id of simulation to load, just write the folder
7 # name here
8
9 # Layer interface positions
10 # defined as decimals
11 # ex: [0.4,0.6,1] will give two boundaries at 0.4 from the bottom
12 # and 0.6 from the bottom. ie 3 layers
13 # include 1 in the end as roof.
14 interfaces = [0,0.4,0.6,1]
15
16 # mechanical properties for each layer
17 youngs_modulus = [2e7,2e7,2e7] # elastic modulus
18 poissons_ratio = [0.185,0.185,0.185] # shear stiffness ratio
19 tensile_strength = [0.3,0.01,0.3] # strength of bonds between grains
20 shear_strength = [0.3,0.01,0.3] # shear stregth of bonds
21 contact_dynamic_friction = [0.4,0.01,0.4] # friction between grains
22 rotating = [true,true,true] # can grains rotate or not
23 color = [1,2,1]
24
25 #carpet_youngs_modulus = 2e7
26 #carpet_poissons_ratio = 0.185
27 #carpet_tensile_strength = 1e16
28 #carpet_contact_dynamic_friction = 0.4
29 #carpet_rotating = true
30 #carpet_shear_strength = 1e16
31
32 sim = Granular.readSimulation("$(id)/comp.jld2")
33 SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
34
35 sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
36
37 Granular.zeroKinematics!(sim) # end any movement
38
39 y_top = -Inf
40 for grain in sim.grains
41 grain.contact_viscosity_normal = 0
42 if y_top < grain.lin_pos[2] + grain.contact_radius
43 global y_top = grain.lin_pos[2] + grain.contact_radius
44 end
45 end
46
47 y_bot = Inf
48 for grain in sim.grains
49 if y_bot > grain.lin_pos[2] - grain.contact_radius
50 global y_bot = grain.lin_pos[2] - grain.contact_radius
51 end
52 end
53
54 # quick fix to make the color = 0 flag for grains belonging to the carpet.
55 # this should be done in the newer versions of init_basin.jl instead
56
57 for grain in sim.grains
58 if grain.lin_pos[2] == -0.05
59 grain.color = 0
60 else
61 grain.color = 1
62 end
63 end
64
65
66 h = y_top-y_bot #depth of basin
67
68 interfaces *= h
69
70 for grain in sim.grains
71
72 for i = 2:size(interfaces,1)
73
74 if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1] && grain.color != 0
75
76 grain.youngs_modulus = youngs_modulus[i-1]
77 grain.poissons_ratio = poissons_ratio[i-1]
78 grain.tensile_strength = tensile_strength[i-1]
79 grain.shear_strength = shear_strength[i-1]
80 grain.contact_dynamic_friction = contact_dynamic_friction[i-1]
81 grain.rotating = rotating[i-1]
82 grain.color = color[i-1]
83
84 end
85 end
86 end
87
88 # Create the bonds between grains by expanding all grains by a small amount
89 # then search and establish contacts and then reduce the size of the grains again
90
91 size_increasing_factor = 1.10 # factor by which contact radius should be increased
92 # to search for contacts
93 increase_array = []
94
95 #increase the contact radius
96 for grain in sim.grains
97 if grain.color != 0
98 contact_radius_increase = (grain.contact_radius*size_increasing_factor)-grain.contact_radius
99 grain.contact_radius += contact_radius_increase
100 append!(increase_array,contact_radius_increase)
101 elseif grain.color == 0
102 append!(increase_array,0)
103 end
104 end
105
106 Granular.findContacts!(sim,method="ocean grid")
107 #Granular.findContactsAllToAll!(sim) # find the grain contacts
108 #Granular.run!(sim,single_step=true)
109
110 #reduce the contact radius again
111 for i = 1:size(sim.grains,1)
112 sim.grains[i].contact_radius -= increase_array[i]
113 end
114
115
116 cd("$id")
117 sim.id = "layered"
118
119 #Granular.resetTime!(sim)
120 #Granular.setTotalTime!(sim,0.5)
121 #Granular.run!(sim)
122
123 Granular.resetTime!(sim)
124 Granular.run!(sim,single_step=true)
125
126 cd("..")
127
128 Granular.writeSimulation(sim,
129 filename = "$(id)/layered.jld2")