commit 7e40836be8f35a54b8259fa307dc952f011bf868
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 4 Aug 2021 10:15:58 +0200
first commit
Diffstat:
| A | LICENSE | | | 15 | +++++++++++++++ |
| A | Makefile | | | 257 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | README | | | 1 | + |
| A | compress.jl | | | 262 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | init-assemblage.jl | | | 219 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | normal-stress.gp | | | 22 | ++++++++++++++++++++++ |
6 files changed, 776 insertions(+), 0 deletions(-)
diff --git a/LICENSE b/LICENSE
@@ -0,0 +1,15 @@
+ISC License
+
+Copyright (c) 2021 Anders Damsgaard <anders@adamsgaard.dk>
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
diff --git a/Makefile b/Makefile
@@ -0,0 +1,257 @@
+JULIA=julia
+INITFILE=rotating01-seed1/rotating01-seed1.681.jld2 compress.jl
+
+FAULTING01 =\
+ faulting01-mu_d-0.8-rotate-seed1-data.tsv \
+ faulting01-mu_d-0.8-norotate-seed1-data.tsv \
+ faulting01-mu_d-0.5-rotate-seed1-data.tsv \
+ faulting01-mu_d-0.5-norotate-seed1-data.tsv \
+ faulting01-mu_d-0.3-rotate-seed1-data.tsv \
+ faulting01-mu_d-0.3-norotate-seed1-data.tsv \
+ faulting01-mu_d-0.1-rotate-seed1-data.tsv \
+ faulting01-mu_d-0.1-norotate-seed1-data.tsv \
+
+FAULTING02 =\
+ faulting02-mu_d-0.8-rotate-seed1-data.tsv \
+ faulting02-mu_d-0.8-norotate-seed1-data.tsv \
+ faulting02-mu_d-0.5-rotate-seed1-data.tsv \
+ faulting02-mu_d-0.5-norotate-seed1-data.tsv \
+ faulting02-mu_d-0.3-rotate-seed1-data.tsv \
+ faulting02-mu_d-0.3-norotate-seed1-data.tsv \
+ faulting02-mu_d-0.1-rotate-seed1-data.tsv \
+ faulting02-mu_d-0.1-norotate-seed1-data.tsv \
+
+FAULTING03 =\
+ faulting03-mu_d-0.8-rotate-seed1-data.tsv \
+ faulting03-mu_d-0.8-norotate-seed1-data.tsv \
+ faulting03-mu_d-0.5-rotate-seed1-data.tsv \
+ faulting03-mu_d-0.5-norotate-seed1-data.tsv \
+ faulting03-mu_d-0.3-rotate-seed1-data.tsv \
+ faulting03-mu_d-0.3-norotate-seed1-data.tsv \
+ faulting03-mu_d-0.1-rotate-seed1-data.tsv \
+ faulting03-mu_d-0.1-norotate-seed1-data.tsv \
+
+FAULTING01_TENSILE_STRENGTH = 0.0
+FAULTING01_SHEAR_STRENGTH = 0.0
+
+FAULTING02_TENSILE_STRENGTH = 400e3
+FAULTING02_SHEAR_STRENGTH = 200e3
+
+FAULTING03_TENSILE_STRENGTH = 40e3
+FAULTING03_SHEAR_STRENGTH = 20e3
+
+all: \
+ faulting01-normal-stress.pdf \
+ faulting02-normal-stress.pdf \
+ faulting03-normal-stress.pdf \
+
+faulting03-normal-stress.pdf: ${FAULTING03} normal-stress.gp
+ gnuplot -e "id=03" normal-stress.gp > faulting03-normal-stress.pdf
+
+faulting02-normal-stress.pdf: ${FAULTING02} normal-stress.gp
+ gnuplot -e "id=02" normal-stress.gp > faulting02-normal-stress.pdf
+
+faulting01-normal-stress.pdf: ${FAULTING01} normal-stress.gp
+ gnuplot -e "id=01" normal-stress.gp > faulting01-normal-stress.pdf
+
+faulting01-mu_d-0.8-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.8-rotate >/dev/null
+
+faulting01-mu_d-0.8-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.8-norotate >/dev/null
+
+faulting01-mu_d-0.5-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.5-rotate >/dev/null
+
+faulting01-mu_d-0.5-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.5-norotate >/dev/null
+
+faulting01-mu_d-0.3-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.3-rotate >/dev/null
+
+faulting01-mu_d-0.3-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.3-norotate >/dev/null
+
+faulting01-mu_d-0.1-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.1-rotate >/dev/null
+
+faulting01-mu_d-0.1-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING01_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING01_SHEAR_STRENGTH} \
+ faulting01-mu_d-0.1-norotate >/dev/null
+
+faulting02-mu_d-0.8-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.8-rotate >/dev/null
+
+faulting02-mu_d-0.8-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.8-norotate >/dev/null
+
+faulting02-mu_d-0.5-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.5-rotate >/dev/null
+
+faulting02-mu_d-0.5-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.5-norotate >/dev/null
+
+faulting02-mu_d-0.3-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.3-rotate >/dev/null
+
+faulting02-mu_d-0.3-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.3-norotate >/dev/null
+
+faulting02-mu_d-0.1-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.1-rotate >/dev/null
+
+faulting02-mu_d-0.1-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING02_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING02_SHEAR_STRENGTH} \
+ faulting02-mu_d-0.1-norotate >/dev/null
+
+faulting03-mu_d-0.8-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.8-rotate >/dev/null
+
+faulting03-mu_d-0.8-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.8 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.8-norotate >/dev/null
+
+faulting03-mu_d-0.5-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.5-rotate >/dev/null
+
+faulting03-mu_d-0.5-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.5 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.5-norotate >/dev/null
+
+faulting03-mu_d-0.3-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.3-rotate >/dev/null
+
+faulting03-mu_d-0.3-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.3 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.3-norotate >/dev/null
+
+faulting03-mu_d-0.1-rotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating true \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.1-rotate >/dev/null
+
+faulting03-mu_d-0.1-norotate-seed1-data.txt: ${INITFILE}
+ ${JULIA} compress.jl \
+ --rotating false \
+ --mu_d 0.1 \
+ --tensile_strength ${FAULTING03_TENSILE_STRENGTH} \
+ --shear_strength ${FAULTING03_SHEAR_STRENGTH} \
+ faulting03-mu_d-0.1-norotate >/dev/null
+
+# transpose(1) from git://src.adamsgaard.dk/mathtools
+.txt.tsv:
+ transpose < $< > $@
+
+${INITFILE}:
+ ${JULIA} init-assemblage.jl --rotating true rotating01
+
+.SUFFIXES: .txt .tsv
+.PHONY: all
diff --git a/README b/README
@@ -0,0 +1 @@
+To run, type 'make'.
diff --git a/compress.jl b/compress.jl
@@ -0,0 +1,262 @@
+#/usr/bin/env julia
+ENV["MPLBACKEND"] = "Agg"
+import Granular
+import PyPlot
+import ArgParse
+using Random
+import DelimitedFiles
+
+let
+render = false
+
+function parse_command_line()
+ s = ArgParse.ArgParseSettings()
+ ArgParse.@add_arg_table s begin
+ "--E"
+ help = "Young's modulus [Pa]"
+ arg_type = Float64
+ default = 2e7
+ "--nu"
+ help = "Poisson's ratio [-]"
+ arg_type = Float64
+ default = 0.285
+ "--mu_d"
+ help = "dynamic friction coefficient [-]"
+ arg_type = Float64
+ default = 0.3
+ "--tensile_strength"
+ help = "the maximum tensile strength [Pa]"
+ arg_type = Float64
+ default = 400e3
+ "--tensile_strength_std_dev"
+ help = "standard deviation of the maximum tensile strength [Pa]"
+ arg_type = Float64
+ default = 0.0
+ "--shear_strength"
+ help = "the maximum shear strength [Pa]"
+ arg_type = Float64
+ default = 200e3
+ "--fracture_toughness"
+ help = "fracture toughness [Pa*m^{1/2}]"
+ arg_type = Float64
+ default = 1e20
+ #default = 1285e3
+ "--r_min"
+ help = "min. grain radius [m]"
+ arg_type = Float64
+ default = 20.0
+ "--r_max"
+ help = "max. grain radius [m]"
+ arg_type = Float64
+ default = 80.0
+ "--thickness"
+ help = "grain thickness [m]"
+ arg_type = Float64
+ default = 1.
+ "--rotating"
+ help = "allow the grains to rotate"
+ arg_type = Bool
+ default = true
+ "--compressive_velocity"
+ help = "compressive velocity [m/s]"
+ arg_type = Float64
+ default = 0.1
+ "--seed"
+ help = "seed for the pseudo RNG"
+ arg_type = Int
+ default = 1
+ "id"
+ help = "simulation id"
+ required = true
+ end
+ return ArgParse.parse_args(s)
+end
+
+function report_args(parsed_args)
+ println("Parsed args:")
+ for (arg,val) in parsed_args
+ println(" $arg => $val")
+ end
+end
+
+function run_simulation(id::String,
+ E::Float64,
+ nu::Float64,
+ mu_d::Float64,
+ tensile_strength::Float64,
+ tensile_strength_std_dev::Float64,
+ shear_strength::Float64,
+ fracture_toughness::Float64,
+ r_min::Float64,
+ r_max::Float64,
+ thickness::Float64,
+ rotating::Bool,
+ compressive_velocity::Float64,
+ seed::Int)
+
+ ############################################################################
+ #### SIMULATION PARAMETERS
+ ############################################################################
+
+ # Grain mechanical properties
+ youngs_modulus = E # Elastic modulus [Pa]
+ poissons_ratio = nu # Shear-stiffness ratio [-]
+ contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
+
+ # Simulation duration [s]
+ t_total = 60.0 * 60.0 * 2.0
+
+ # Temporal interval between output files
+ file_dt = 18.0
+
+ ############################################################################
+ #### Read prior state
+ ############################################################################
+ sim = Granular.readSimulation("rotating01-seed1/rotating01-seed1.681.jld2")
+ sim.id = id
+ Granular.resetTime!(sim)
+ Granular.zeroKinematics!(sim)
+ Granular.removeSimulationFiles(sim)
+
+ # Set grain mechanical properties
+ for grain in sim.grains
+ grain.youngs_modulus = youngs_modulus
+ grain.poissons_ratio = poissons_ratio
+ grain.tensile_strength = abs(tensile_strength +
+ randn()*tensile_strength_std_dev)
+ grain.shear_strength = abs(shear_strength +
+ randn()*tensile_strength_std_dev)
+ grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
+ grain.contact_dynamic_friction = contact_dynamic_friction
+ grain.fracture_toughness = fracture_toughness
+ grain.rotating = rotating
+ grain.ocean_drag_coeff_horiz = 1.6e-4
+ grain.ocean_drag_coeff_vert = 0.14
+ end
+
+ x_max = -Inf
+ x_min = +Inf
+ y_max = -Inf
+ y_min = +Inf
+ r_max = 0.0
+ for grain in sim.grains
+ r = grain.contact_radius
+
+ if x_min > grain.lin_pos[1] - r
+ x_min = grain.lin_pos[1] - r
+ end
+
+ if x_max < grain.lin_pos[1] + r
+ x_max = grain.lin_pos[1] + r
+ end
+
+ if y_min > grain.lin_pos[2] - r
+ y_min = grain.lin_pos[2] - r
+ end
+
+ if y_max < grain.lin_pos[2] + r
+ y_max = grain.lin_pos[2] + r
+ end
+
+ if r_max < r
+ r_max = r
+ end
+ end
+
+ dx = r_max * 2.0
+ L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx]
+ n = convert(Vector{Int}, floor.(L./dx))
+ sim.ocean = Granular.createRegularOceanGrid(n, L,
+ origo=[x_min - L[1]/3.0, y_min],
+ time=[0.], name="resized")
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
+
+ sim.walls[1].bc = "velocity"
+ sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total
+
+ # fix grains adjacent to walls to create no-slip BC
+ fixed_thickness = 2.0 * r_max
+ for grain in sim.grains
+ if grain.lin_pos[2] <= y_min + fixed_thickness
+ grain.fixed = true
+ grain.color = 1
+ elseif grain.lin_pos[2] >= y_max - fixed_thickness
+ grain.allow_y_acc = true
+ grain.fixed = true
+ grain.color = 1
+ end
+ end
+
+ Granular.setTimeStep!(sim)
+ Granular.setTotalTime!(sim, t_total)
+ Granular.setOutputFileInterval!(sim, file_dt)
+ render && Granular.plotGrains(sim, verbose=true)
+
+ ############################################################################
+ #### RUN THE SIMULATION
+ ############################################################################
+
+ # Run the consolidation experiment, and monitor top wall position over
+ # time
+ time = Float64[]
+ compaction = Float64[]
+ effective_normal_stress = Float64[]
+ while sim.time < sim.time_total
+
+ for i=1:100
+ 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
+ Granular.writeSimulation(sim)
+
+ defined_normal_stress = ones(length(effective_normal_stress)) *
+ 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 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()
+ DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction,
+ effective_normal_stress,
+ defined_normal_stress])
+
+ render && Granular.render(sim, trim=false, animation=true)
+end
+
+function main()
+ parsed_args = parse_command_line()
+ report_args(parsed_args)
+
+ seed = parsed_args["seed"]
+ run_simulation(parsed_args["id"] * "-seed" * string(seed),
+ parsed_args["E"],
+ parsed_args["nu"],
+ parsed_args["mu_d"],
+ parsed_args["tensile_strength"],
+ parsed_args["tensile_strength_std_dev"],
+ parsed_args["shear_strength"],
+ parsed_args["fracture_toughness"],
+ parsed_args["r_min"],
+ parsed_args["r_max"],
+ parsed_args["thickness"],
+ parsed_args["rotating"],
+ parsed_args["compressive_velocity"],
+ seed)
+end
+
+main()
+end
diff --git a/init-assemblage.jl b/init-assemblage.jl
@@ -0,0 +1,219 @@
+#/usr/bin/env julia
+ENV["MPLBACKEND"] = "Agg"
+import Granular
+import PyPlot
+import ArgParse
+using Random
+import DelimitedFiles
+
+let
+render = false
+
+function parse_command_line()
+ s = ArgParse.ArgParseSettings()
+ ArgParse.@add_arg_table s begin
+ "--E"
+ help = "Young's modulus [Pa]"
+ arg_type = Float64
+ default = 2e7
+ "--nu"
+ help = "Poisson's ratio [-]"
+ arg_type = Float64
+ default = 0.285
+ "--mu_d"
+ help = "dynamic friction coefficient [-]"
+ arg_type = Float64
+ default = 0.3
+ "--tensile_strength"
+ help = "the maximum tensile strength [Pa]"
+ arg_type = Float64
+ default = 400e3
+ "--tensile_strength_std_dev"
+ help = "standard deviation of the maximum tensile strength [Pa]"
+ arg_type = Float64
+ default = 0.0
+ "--shear_strength"
+ help = "the maximum shear strength [Pa]"
+ arg_type = Float64
+ default = 200e3
+ "--fracture_toughness"
+ help = "fracture toughness [Pa*m^{1/2}]"
+ arg_type = Float64
+ default = 1e20
+ #default = 1285e3
+ "--r_min"
+ help = "min. grain radius [m]"
+ arg_type = Float64
+ default = 20.0
+ "--r_max"
+ help = "max. grain radius [m]"
+ arg_type = Float64
+ default = 80.0
+ "--thickness"
+ help = "grain thickness [m]"
+ arg_type = Float64
+ default = 1.
+ "--rotating"
+ help = "allow the grains to rotate"
+ arg_type = Bool
+ default = true
+ "--compressive_velocity"
+ help = "compressive velocity [m/s]"
+ arg_type = Float64
+ default = 0.1
+ "--seed"
+ help = "seed for the pseudo RNG"
+ arg_type = Int
+ default = 1
+ "id"
+ help = "simulation id"
+ required = true
+ end
+ return ArgParse.parse_args(s)
+end
+
+function report_args(parsed_args)
+ println("Parsed args:")
+ for (arg,val) in parsed_args
+ println(" $arg => $val")
+ end
+end
+
+function run_simulation(id::String,
+ E::Float64,
+ nu::Float64,
+ mu_d::Float64,
+ tensile_strength::Float64,
+ tensile_strength_std_dev::Float64,
+ shear_strength::Float64,
+ fracture_toughness::Float64,
+ r_min::Float64,
+ r_max::Float64,
+ thickness::Float64,
+ rotating::Bool,
+ compressive_velocity::Float64,
+ seed::Int)
+
+ ############################################################################
+ #### SIMULATION PARAMETERS
+ ############################################################################
+
+ # Common simulation identifier
+ id_prefix = id
+
+ # Grain mechanical properties
+ youngs_modulus = E # Elastic modulus [Pa]
+ poissons_ratio = nu # Shear-stiffness ratio [-]
+ contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
+
+ # Simulation duration [s]
+ t_total = 60.0 * 60.0 * 4.0
+
+ # Temporal interval between output files
+ file_dt = 18.0
+
+ # Misc parameters
+ water_density = 1000.
+
+ # World size [m]
+ L = [5e3, 2e4, 10.0]
+
+ Random.seed!(seed)
+
+ ############################################################################
+ #### Create ice floes
+ ############################################################################
+ sim = Granular.createSimulation("$(id_prefix)")
+ Granular.removeSimulationFiles(sim)
+
+ # Create box edges of ice floes with r_min radius, moving inward
+ r_walls = r_min
+
+ # Create a grid for contact searching spanning the extent of the grains
+ n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2]
+ sim.ocean = Granular.createRegularOceanGrid(n, L)
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
+ Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
+
+ # Add unconstrained ice floes
+ Granular.irregularPacking!(sim,
+ radius_max=r_max, radius_min=r_min,
+ thickness=thickness,
+ #binary_radius_search=true,
+ seed=seed)
+
+ # Set grain mechanical properties
+ for grain in sim.grains
+ grain.youngs_modulus = youngs_modulus
+ grain.poissons_ratio = poissons_ratio
+ grain.tensile_strength = abs(tensile_strength +
+ randn()*tensile_strength_std_dev)
+ grain.shear_strength = abs(shear_strength +
+ randn()*tensile_strength_std_dev)
+ grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
+ grain.contact_dynamic_friction = contact_dynamic_friction
+ grain.fracture_toughness = fracture_toughness
+ grain.rotating = rotating
+ grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
+ end
+
+ # Create GSD plot to signify that simulation is complete
+ Granular.writeSimulation(sim)
+ render && Granular.plotGrainSizeDistribution(sim)
+
+ # Add a dynamic wall to the top which adds a normal stress downwards.
+ # The normal of this wall is downwards, and we place it at the top of
+ # the granular assemblage. Here, the fluid drag is also disabled
+ y_top = -Inf
+ for grain in sim.grains
+ if y_top < grain.lin_pos[2] + grain.contact_radius
+ y_top = grain.lin_pos[2] + grain.contact_radius
+ end
+ end
+ Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
+ bc="normal stress",
+ normal_stress=-20e3)
+ sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass
+ @info("Placing top wall at y=$y_top")
+
+ # Resize the grid to span the current state
+ Granular.fitGridToGrains!(sim, sim.ocean)
+
+ Granular.setTimeStep!(sim)
+ Granular.setTotalTime!(sim, t_total)
+ Granular.setOutputFileInterval!(sim, file_dt)
+ render && Granular.plotGrains(sim, verbose=true)
+
+ ############################################################################
+ #### RUN THE SIMULATION
+ ############################################################################
+
+ Granular.run!(sim)
+ Granular.writeSimulation(sim)
+
+ render && Granular.render(sim, trim=false, animation=true)
+end
+
+function main()
+ parsed_args = parse_command_line()
+ report_args(parsed_args)
+
+ seed = parsed_args["seed"]
+ run_simulation(parsed_args["id"] * "-seed" * string(seed),
+ parsed_args["E"],
+ parsed_args["nu"],
+ parsed_args["mu_d"],
+ parsed_args["tensile_strength"],
+ parsed_args["tensile_strength_std_dev"],
+ parsed_args["shear_strength"],
+ parsed_args["fracture_toughness"],
+ parsed_args["r_min"],
+ parsed_args["r_max"],
+ parsed_args["thickness"],
+ parsed_args["rotating"],
+ parsed_args["compressive_velocity"],
+ seed)
+end
+
+main()
+end
diff --git a/normal-stress.gp b/normal-stress.gp
@@ -0,0 +1,22 @@
+#!/usr/bin/env gnuplot
+
+reset
+
+set terminal pdfcairo color size 17.8 cm, 9.45 cm font ",10"
+
+set multiplot layout 2,1 #\
+ #margin 0.07,0.92,0.25,0.86 \
+ #spacing 0.17,0.0
+
+set key bottom right #samplen 0.9
+
+set xlabel "Compressive strain {/:Italic e} [-]"
+set ylabel "Normal stress {/:Italic N} [kPa]"
+
+set yrange [0:60]
+
+file(mu, rotate) = sprintf("faulting%02d-mu_d-%s-%s-seed1-data.tsv", id, mu, rotate)
+label(mu, rotate) = sprintf("{/symbol m}_d = %s, %s", mu, rotate)
+
+plot for [mu in "0.1 0.3 0.5 0.8"] file(mu, "rotate") u ($0 / 4246 * 0.2):($3 / 1000) w l t label(mu, "rotating"),
+ plot for [mu in "0.1 0.3 0.5 0.8"] file(mu, "norotate") u ($0 / 4246 * 0.2):($3 / 1000) w l t label(mu, "not rotating")
