seaice-experiments

sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments # fast
git clone https://src.adamsgaard.dk/seaice-experiments.git # slow
Log | Files | Refs | README | LICENSE Back to index

ridging_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")