Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit e6423e746bc3d5ee39a065427f18b4be7a8021b3
parent 6682dda64ebfab3793c6d916f258f5925892df6f
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Thu, 11 May 2017 16:45:39 -0400

add angular drag from ocean vorticity

Diffstat:
Msrc/grid.jl | 39++++++++++++++++++++++++++++++++++++++-
Msrc/ocean.jl | 27+++++++++++++++++++++++++--
Mtest/grid.jl | 15+++++++++++++++
Mtest/ocean.jl | 142++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------
4 files changed, 201 insertions(+), 22 deletions(-)

diff --git a/src/grid.jl b/src/grid.jl @@ -7,7 +7,8 @@ south-west (-x, -y)-facing corner. # Arguments * `field::Array{Float64, 4}`: a scalar field to interpolate from -* `p::float`: point position +* `x_tilde::float`: x point position [0;1] +* `y_tilde::float`: y point position [0;1] * `i::Int`: i-index of cell containing point * `j::Int`: j-index of scalar field to interpolate from * `it::Int`: time step from scalar field to interpolate from @@ -30,6 +31,42 @@ function bilinearInterpolation(field::Array{Float64, 4}, field[i, j, k, it]*(1. - x_tilde))*(1. - y_tilde) end +""" + curl(ocean, x_tilde, y_tilde, i, j, k, it) + +Use bilinear interpolation to interpolate curl value for a staggered velocity +grid to an arbitrary position in a cell. Assumes south-west convention, i.e. +(i,j) is located at the south-west (-x, -y)-facing corner. + +# Arguments +* `ocean::Ocean`: grid for which to determine curl +* `x_tilde::float`: x point position [0;1] +* `y_tilde::float`: y point position [0;1] +* `i::Int`: i-index of cell containing point +* `j::Int`: j-index of scalar field to interpolate from +* `it::Int`: time step from scalar field to interpolate from +""" +function curl(ocean::Ocean, + x_tilde::Float64, + y_tilde::Float64, + i::Int, + j::Int, + k::Int, + it::Int) + + sw, se, ne, nw = getCellCornerCoordinates(ocean, i, j) + sw_se = norm(sw - se) + se_ne = norm(se - ne) + nw_ne = norm(nw - ne) + sw_nw = norm(sw - nw) + + return ( + ((ocean.v[i+1, j , k,it] - ocean.v[i , j , k,it])/sw_se*(1. - y_tilde) + + ((ocean.v[i+1, j+1, k,it] - ocean.v[i , j+1, k,it])/nw_ne)*y_tilde) - + ((ocean.u[i , j+1, k,it] - ocean.u[i , j , k,it])/sw_nw*(1. - x_tilde) + + ((ocean.u[i+1, j+1, k,it] - ocean.u[i+1, j , k,it])/se_ne)*x_tilde)) +end + export sortIceFloesInOceanGrid! """ Find ice-floe positions in ocean grid, based on their center positions. diff --git a/src/ocean.jl b/src/ocean.jl @@ -245,7 +245,8 @@ end export addOceanDrag! """ -Add Stokes-type drag from velocity difference between ocean and all ice floes. +Add drag from linear and angular velocity difference between ocean and all ice +floes. """ function addOceanDrag!(simulation::Simulation) if typeof(simulation.ocean.input_file) == Bool @@ -276,8 +277,10 @@ function addOceanDrag!(simulation::Simulation) u_local = bilinearInterpolation(u, x_tilde, y_tilde, i, j, k, 1) v_local = bilinearInterpolation(v, x_tilde, y_tilde, i, j, k, 1) + vel_curl = curl(simulation.ocean, x_tilde, y_tilde, i, j, k, 1) applyOceanDragToIceFloe!(ice_floe, u_local, v_local) + applyOceanVorticityToIceFloe!(ice_floe, vel_curl) end end @@ -297,6 +300,26 @@ function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical, width = ice_floe.areal_radius*2. ice_floe.force += - rho_o * (.5*c_o_v*width*draft*freeboard + c_o_h*length*width) * + rho_o * (.5*c_o_v*width*draft + c_o_h*length*width) * ([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel) end + +export applyOceanVorticityToIceFloe! +""" +Add Stokes-type torque from angular velocity difference between ocean and a +single ice floe. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama +and Boucher, 1999. +""" +function applyOceanVorticityToIceFloe!(ice_floe::IceFloeCylindrical, + ocean_curl::float) + 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)* + abs(.5*ocean_curl - ice_floe.ang_vel)*(.5*ocean_curl - ice_floe.ang_vel) +end diff --git a/test/grid.jl b/test/grid.jl @@ -137,3 +137,18 @@ SeaIce.addOceanDrag!(sim) @test sim.ice_floes[1].force[2] > 0. @test sim.ice_floes[2].force[1] < 0. @test sim.ice_floes[2].force[2] > 0. + +info("Testing curl function") +ocean.u[1, 1, 1, 1] = 1.0 +ocean.u[2, 1, 1, 1] = 1.0 +ocean.u[2, 2, 1, 1] = 0.0 +ocean.u[1, 2, 1, 1] = 0.0 +ocean.v[:, :, 1, 1] = 0.0 +@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1) > 0. + +ocean.u[1, 1, 1, 1] = 0.0 +ocean.u[2, 1, 1, 1] = 0.0 +ocean.u[2, 2, 1, 1] = 1.0 +ocean.u[1, 2, 1, 1] = 1.0 +ocean.v[:, :, 1, 1] = 0.0 +@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1) < 0. diff --git a/test/ocean.jl b/test/ocean.jl @@ -6,22 +6,126 @@ info("#### $(basename(@__FILE__)) ####") info("Testing regular grid generation") -ocean = SeaIce.createRegularOceanGrid([6, 6, 6], [1., 1., 1.]) -@test size(ocean.xq) == (7, 7) -@test size(ocean.yq) == (7, 7) -@test size(ocean.xh) == (6, 6) -@test size(ocean.yh) == (6, 6) -@test ocean.xq[1, :, 1] ≈ zeros(7) -@test ocean.xq[4, :, 1] ≈ .5*ones(7) -@test ocean.xq[end, :, 1] ≈ 1.*ones(7) -@test ocean.yq[:, 1, 1] ≈ zeros(7) -@test ocean.yq[:, 4, 1] ≈ .5*ones(7) -@test ocean.yq[:, end, 1] ≈ 1.*ones(7) -@test size(ocean.u) == (7, 7, 6, 1) -@test size(ocean.v) == (7, 7, 6, 1) -@test size(ocean.h) == (7, 7, 6, 1) -@test size(ocean.e) == (7, 7, 6, 1) -@test ocean.u ≈ zeros(7, 7, 6, 1) -@test ocean.v ≈ zeros(7, 7, 6, 1) -@test ocean.h ≈ zeros(7, 7, 6, 1) -@test ocean.e ≈ zeros(7, 7, 6, 1) +sim = SeaIce.createSimulation() +sim.ocean = SeaIce.createRegularOceanGrid([6, 6, 6], [1., 1., 1.]) +@test size(sim.ocean.xq) == (7, 7) +@test size(sim.ocean.yq) == (7, 7) +@test size(sim.ocean.xh) == (6, 6) +@test size(sim.ocean.yh) == (6, 6) +@test sim.ocean.xq[1, :, 1] ≈ zeros(7) +@test sim.ocean.xq[4, :, 1] ≈ .5*ones(7) +@test sim.ocean.xq[end, :, 1] ≈ 1.*ones(7) +@test sim.ocean.yq[:, 1, 1] ≈ zeros(7) +@test sim.ocean.yq[:, 4, 1] ≈ .5*ones(7) +@test sim.ocean.yq[:, end, 1] ≈ 1.*ones(7) +@test size(sim.ocean.u) == (7, 7, 6, 1) +@test size(sim.ocean.v) == (7, 7, 6, 1) +@test size(sim.ocean.h) == (7, 7, 6, 1) +@test size(sim.ocean.e) == (7, 7, 6, 1) +@test sim.ocean.u ≈ zeros(7, 7, 6, 1) +@test sim.ocean.v ≈ zeros(7, 7, 6, 1) +@test sim.ocean.h ≈ zeros(7, 7, 6, 1) +@test sim.ocean.e ≈ zeros(7, 7, 6, 1) + +info("Testing velocity drag interaction (static ocean)") +SeaIce.addIceFloeCylindrical(sim, [.5, .5], .25, .1) +SeaIce.setTotalTime!(sim, 5.) +SeaIce.setTimeStep!(sim) +sim_init = deepcopy(sim) +sim.ice_floes[1].lin_vel[1] = 0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after +@test E_kin_lin_init > E_kin_lin_final # linear velocity lost due to ocean drag + +info("Testing velocity drag interaction (static ice floe)") +sim = deepcopy(sim_init) +sim.ocean.v[:, :, 1, 1] = 0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after +@test E_kin_lin_init < E_kin_lin_final # linear vel. gained due to ocean drag + +info("Testing vortex interaction between (static ocean)") +sim = deepcopy(sim_init) +sim.ice_floes[1].ang_vel = 0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test E_kin_rot_init > E_kin_rot_final # energy lost to ocean +@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation +@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation +@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained + +info("Testing vortex interaction between (static ice floe)") +sim = deepcopy(sim_init) +sim.ocean = SeaIce.createRegularOceanGrid([1, 1, 1], [1., 1., 1.]) +sim.ice_floes[1].lin_pos[1] = 0.5 +sim.ice_floes[1].lin_pos[2] = 0.5 +sim.ocean.v[1, :, 1, 1] = -0.1 +sim.ocean.v[2, :, 1, 1] = 0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation +@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation +@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex +@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained + +sim = deepcopy(sim_init) +sim.ocean = SeaIce.createRegularOceanGrid([1, 1, 1], [1., 1., 1.]) +sim.ice_floes[1].lin_pos[1] = 0.5 +sim.ice_floes[1].lin_pos[2] = 0.5 +sim.ocean.v[1, :, 1, 1] = 0.1 +sim.ocean.v[2, :, 1, 1] = -0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test sim.ice_floes[1].ang_vel < 0. # check angular velocity orientation +@test sim.ice_floes[1].ang_pos < 0. # check angular position orientation +@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex +@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained + +sim = deepcopy(sim_init) +sim.ocean = SeaIce.createRegularOceanGrid([1, 1, 1], [1., 1., 1.]) +sim.ice_floes[1].lin_pos[1] = 0.5 +sim.ice_floes[1].lin_pos[2] = 0.5 +sim.ocean.u[:, 1, 1, 1] = -0.1 +sim.ocean.u[:, 2, 1, 1] = 0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test sim.ice_floes[1].ang_vel < 0. # check angular velocity orientation +@test sim.ice_floes[1].ang_pos < 0. # check angular position orientation +@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex +@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained + +sim = deepcopy(sim_init) +sim.ocean = SeaIce.createRegularOceanGrid([1, 1, 1], [1., 1., 1.]) +sim.ice_floes[1].lin_pos[1] = 0.5 +sim.ice_floes[1].lin_pos[2] = 0.5 +sim.ocean.u[:, 1, 1, 1] = 0.1 +sim.ocean.u[:, 2, 1, 1] = -0.1 +E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +SeaIce.run!(sim, verbose=false) +E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim) +E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim) +@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation +@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation +@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex +@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained