commit b4856096f45707589baa494ad735de2128f7ec35
parent d9942595aec41f37d766c59d041229d6f738e0c0
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed, 24 May 2017 11:54:04 -0400
save additional contact parameters and drag coefficients in ice floe type
Diffstat:
5 files changed, 102 insertions(+), 11 deletions(-)
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -44,6 +44,17 @@ type IceFloeCylindrical
     contact_static_friction::float
     contact_dynamic_friction::float
 
+    youngs_modulus::float
+    poissons_ratio::float
+    tensile_strength::float
+    compressive_strength_prefactor::float
+
+    # Ocean/atmosphere interaction parameters
+    ocean_drag_coeff_vert::float
+    ocean_drag_coeff_horiz::float
+    atmos_drag_coeff_vert::float
+    atmos_drag_coeff_horiz::float
+
     # Interaction
     pressure::float
     n_contacts::Int
@@ -94,6 +105,16 @@ type IceFloeArrays
     contact_static_friction
     contact_dynamic_friction
 
+    youngs_modulus
+    poissons_ratio
+    tensile_strength
+    compressive_strength_prefactor
+
+    ocean_drag_coeff_vert
+    ocean_drag_coeff_horiz
+    atmos_drag_coeff_vert
+    atmos_drag_coeff_horiz
+
     pressure
     n_contacts
 end
diff --git a/src/icefloe.jl b/src/icefloe.jl
@@ -27,6 +27,15 @@ function addIceFloeCylindrical(simulation::Simulation,
                                contact_viscosity_tangential::float = 0.,
                                contact_static_friction::float = 0.4,
                                contact_dynamic_friction::float = 0.4,
+                               youngs_modulus::float = 2e9,  # Hopkins 2004
+                               poissons_ratio::float = 0.185,  # Hopkins 2004
+                               tensile_strength::float = 500e3,  # Hopkins 2004
+                               compressive_strength_prefactor::float = 1285e3,  
+                                   # Hopkins 2004
+                               ocean_drag_coeff_vert::float = 0.85, # H&C 2011
+                               ocean_drag_coeff_horiz::float = 5e-4, # H&C 2011
+                               atmos_drag_coeff_vert::float = 0.4, # H&C 2011
+                               atmos_drag_coeff_horiz::float = 2.5e-4, # H&C2011
                                pressure::float = 0.,
                                fixed::Bool = false,
                                rotating::Bool = true,
@@ -102,6 +111,16 @@ function addIceFloeCylindrical(simulation::Simulation,
                                  contact_static_friction,
                                  contact_dynamic_friction,
 
+                                 youngs_modulus,
+                                 poissons_ratio,
+                                 tensile_strength,
+                                 compressive_strength_prefactor,
+
+                                 ocean_drag_coeff_vert,
+                                 ocean_drag_coeff_horiz,
+                                 atmos_drag_coeff_vert,
+                                 atmos_drag_coeff_horiz,
+
                                  pressure,
                                  n_contacts,
                                  ocean_grid_pos,
@@ -199,6 +218,16 @@ function convertIceFloeDataToArrays(simulation::Simulation)
                           Array(Float64, length(simulation.ice_floes)),
 
                           Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+
+                          Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+                          Array(Float64, length(simulation.ice_floes)),
+
+                          Array(Float64, length(simulation.ice_floes)),
                           Array(Int, length(simulation.ice_floes))
                          )
 
@@ -244,6 +273,21 @@ function convertIceFloeDataToArrays(simulation::Simulation)
         ifarr.contact_dynamic_friction[i] = 
             simulation.ice_floes[i].contact_dynamic_friction
 
+        ifarr.youngs_modulus[i] = simulation.ice_floes[i].youngs_modulus
+        ifarr.poissons_ratio[i] = simulation.ice_floes[i].poissons_ratio
+        ifarr.tensile_strength[i] = simulation.ice_floes[i].tensile_strength
+        ifarr.compressive_strength_prefactor[i] = 
+            simulation.ice_floes[i].compressive_strength_prefactor
+
+        ifarr.ocean_drag_coeff_vert[i] = 
+            simulation.ice_floes[i].ocean_drag_coeff_vert
+        ifarr.ocean_drag_coeff_horiz[i] = 
+            simulation.ice_floes[i].ocean_drag_coeff_horiz
+        ifarr.atmos_drag_coeff_vert[i] = 
+            simulation.ice_floes[i].atmos_drag_coeff_vert
+        ifarr.atmos_drag_coeff_horiz[i] = 
+            simulation.ice_floes[i].atmos_drag_coeff_horiz
+
         ifarr.pressure[i] = simulation.ice_floes[i].pressure
 
         ifarr.n_contacts[i] = simulation.ice_floes[i].n_contacts
@@ -286,10 +330,20 @@ function printIceFloeInfo(f::IceFloeCylindrical)
 
     println("  k_n:     $(f.contact_stiffness_normal) N/m")
     println("  k_t:     $(f.contact_stiffness_tangential) N/m")
-    println("  gamma_n: $(f.contact_viscosity_normal) N/(m/s)")
-    println("  gamma_t: $(f.contact_viscosity_tangential) N/(m/s)")
-    println("  mu_s:    $(f.contact_static_friction)")
-    println("  mu_d:    $(f.contact_dynamic_friction)\n")
+    println("  γ_n: $(f.contact_viscosity_normal) N/(m/s)")
+    println("  γ_t: $(f.contact_viscosity_tangential) N/(m/s)")
+    println("  μ_s:    $(f.contact_static_friction)")
+    println("  μ_d:    $(f.contact_dynamic_friction)\n")
+
+    println("  E:      $(f.youngs_modulus) Pa")
+    println("  ν:      $(f.poissons_ratio)")
+    println("  σ_t:    $(f.tensile_strength) Pa")
+    println("  c(σ_c): $(f.compressive_strength_prefactor) m^0.5 Pa\n")
+
+    println("  c_o_v:  $(f.ocean_drag_coeff_vert)")
+    println("  c_o_h:  $(f.ocean_drag_coeff_horiz)")
+    println("  c_a_v:  $(f.atmos_drag_coeff_vert)")
+    println("  c_a_h:  $(f.atmos_drag_coeff_horiz)\n")
 
     println("  pressure:   $(f.pressure) Pa")
     println("  n_contacts: $(f.n_contacts)")
diff --git a/src/io.jl b/src/io.jl
@@ -93,6 +93,24 @@ function writeIceFloeVTK(simulation::Simulation,
     WriteVTK.vtk_point_data(vtkfile, ifarr.contact_dynamic_friction,
                             "Contact friction (dynamic) [-]")
 
+    WriteVTK.vtk_point_data(vtkfile, ifarr.youngs_modulus,
+                            "Young's modulus [Pa]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.poissons_ratio,
+                            "Poisson's ratio [-]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.tensile_strength,
+                            "Tensile strength [Pa]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.compressive_strength_prefactor,
+                            "Compressive strength prefactor [m^0.5 Pa]")
+
+    WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_drag_coeff_vert,
+                            "Ocean drag coefficient (vertical) [-]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.ocean_drag_coeff_horiz,
+                            "Ocean drag coefficient (horizontal) [-]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.atmos_drag_coeff_vert,
+                            "Atmosphere drag coefficient (vertical) [-]")
+    WriteVTK.vtk_point_data(vtkfile, ifarr.atmos_drag_coeff_horiz,
+                            "Atmosphere drag coefficient (horizontal) [-]")
+
     WriteVTK.vtk_point_data(vtkfile, ifarr.pressure,
                             "Contact pressure [Pa]")
 
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -294,13 +294,12 @@ function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical,
     freeboard = .1*ice_floe.thickness  # height above water
     rho_o = 1000.   # ocean density
     draft = ice_floe.thickness - freeboard  # height of submerged thickness
-    c_o_v = .85  # ocean drag coefficient, vertical, Hunke and Comeau 2011
-    c_o_h = 5e-4  # ocean drag coefficient, horizontal
     length = ice_floe.areal_radius*2.
     width = ice_floe.areal_radius*2.
 
     ice_floe.force +=
-        rho_o * (.5*c_o_v*width*draft + c_o_h*length*width) *
+        rho_o * (.5*ice_floe.ocean_drag_coeff_vert*width*draft + 
+        ice_floe.ocean_drag_coeff_horiz*length*width) *
         ([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel)
 end
 
@@ -315,11 +314,10 @@ function applyOceanVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
     freeboard = .1*ice_floe.thickness  # height above water
     rho_o = 1000.   # ocean density
     draft = ice_floe.thickness - freeboard  # height of submerged thickness
-    c_o_v = .85  # ocean drag coefficient, vertical, Hunke and Comeau 2011
-    c_o_h = 5e-4  # ocean drag coefficient, horizontal
 
     ice_floe.torque +=
         pi*ice_floe.areal_radius^4.*rho_o*
-        (ice_floe.areal_radius/5.*c_o_h + draft*c_o_h)*
+        (ice_floe.areal_radius/5.*ice_floe.ocean_drag_coeff_horiz + 
+        draft*ice_floe.ocean_drag_coeff_vert)*
         abs(.5*ocean_curl - ice_floe.ang_vel)*(.5*ocean_curl - ice_floe.ang_vel)
 end
diff --git a/test/vtk.jl b/test/vtk.jl
@@ -26,7 +26,7 @@ else
 end
 
 icefloechecksum = 
-"88daceb1b99c519154b1acdcf8f340967794c552c74ea70c4af8954d8af5296a  " *
+"4885bc7c0eccc5e43a77d370295841fae710465c8bcc120a13681a0947ffbd53  " *
 "test.icefloes.1.vtu\n"
 
 oceanchecksum =