tridging_parameterization2.jl - seaice-experiments - sea ice experiments using Granular.jl
 (HTM) git clone git://src.adamsgaard.dk/seaice-experiments
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tridging_parameterization2.jl (3366B)
       ---
            1 #/usr/bin/env julia
            2 ENV["MPLBACKEND"] = "Agg"
            3 import PyPlot
            4 import LsqFit
            5 import Granular
            6 import DelimitedFiles
            7 import Statistics
            8 
            9 
           10 d = 0.02   # ice particle thickness [m]
           11 L = 200*d  # ice-floe length
           12 E = 3.5e6    # Young's modulus
           13 ν = 0.285  # Poisson's ratio [-]
           14 ρ_w = 1e3  # Water density [kg/m^3]
           15 g = 9.8    # Gravitational acceleration [m/s^2]
           16 
           17 
           18 ################################################################################
           19 #### Compare ridging compressive stress with Granular.jl parameterization
           20 
           21 h1 = Float64[]         # thickness of ice floe 1 [m]
           22 h2 = Float64[]         # thickness of ice floe 2 [m]
           23 comp_vel = Float64[]   # compressive velocity [m/s]
           24 N_max = Float64[]      # compressive stress at yield point [Pa]
           25 τ_max = Float64[]      # shear stress at yield point [Pa]
           26 t_yield = Float64[]    # time of yield failure [t]
           27 d_yield = Float64[]    # compressive distance at yield failure [m]
           28 N_late_avg = Float64[] # averaged post-failure compressive stress [Pa]
           29 τ_late_avg = Float64[] # averaged post-failure shear stress [Pa]
           30 tensile_strength = Float64[] # tensile strength [Pa]
           31 
           32 nx1 = 100
           33 nx2 = 100
           34 ny1 = 10
           35 ny2 = 15
           36 h1 = ny1*d*sin(π/3) # correct thickness for triangular packing
           37 h2 = ny2*d*sin(π/3) # correct thickness for triangular packing
           38 cv = 0.2
           39 sim = Granular.createSimulation()
           40 fracture_toughness = 1.9e6
           41 r1 = nx1*d
           42 r2 = nx2*d
           43 coulomb_friction = 0.4
           44 Granular.addGrainCylindrical!(sim, [0., 0.], r1, h1,
           45                               fracture_toughness=fracture_toughness,
           46                               youngs_modulus=E,
           47                               contact_static_friction=coulomb_friction,
           48                               contact_dynamic_friction=coulomb_friction,
           49                               lin_vel=[cv, 0.],
           50                               fixed=true)
           51 Granular.addGrainCylindrical!(sim, [r1 + r2 + 12*d, 0.0], r2, h2,
           52                               fracture_toughness=fracture_toughness,
           53                               youngs_modulus=E,
           54                               contact_static_friction=coulomb_friction,
           55                               contact_dynamic_friction=coulomb_friction,
           56                               fixed=true)
           57 Granular.setTimeStep!(sim)
           58 compressive_distance = 3.0
           59 Granular.setTotalTime!(sim, compressive_distance/cv)
           60 Granular.removeSimulationFiles(sim)
           61 Granular.setOutputFileInterval!(sim, compressive_distance/cv * (1.0/10.0))
           62 dist = Float64[]
           63 f_n_parameterization = Float64[]
           64 N_parameterization = Float64[]
           65 #println("Peak stress by fit: $(our_strength_model(min(h1[1], h2[1]), fit.param[1])) Pa")
           66 r12 = 2.0*r1*r2/(r1 + r2)
           67 Lz_12 = min(h1, h2)
           68 while sim.time < sim.time_total
           69     for i=1:1
           70         Granular.run!(sim, single_step=true)
           71     end
           72     append!(dist, cv*sim.time)
           73     append!(f_n_parameterization, -sim.grains[1].force[1])
           74     append!(N_parameterization, -sim.grains[1].force[1]/(r12*Lz_12))
           75 end
           76 
           77 PyPlot.figure(figsize=(4,3))
           78 A_exp = ny1*d * d
           79 # Subtract drag against water from data set
           80 #PyPlot.plot(data[1,:].*cv, (data[2,:] - mean(data[2,1:100]))./1e3, "C1-",
           81 #            label="Two-floe experiment")
           82 PyPlot.plot(dist, N_parameterization./1e3, "C2-", label="Plan-view parameterization")
           83 PyPlot.xlabel("Compressive distance [m]")
           84 PyPlot.ylabel("Compressive stress [kPa]")
           85 #PyPlot.legend()
           86 PyPlot.tight_layout()
           87 PyPlot.savefig("ridging_parameterization2.png")
           88 PyPlot.savefig("ridging_parameterization2.pdf")