commit 6a2b205579d371f3ce2760cc445517592956600c
parent b640f26eebec4330354379ec734a996089ec7def
Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 21 Feb 2018 11:28:03 -0500
Keep track of thermal energy on the particle level
Diffstat:
6 files changed, 67 insertions(+), 8 deletions(-)
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -70,6 +70,8 @@ mutable struct GrainCylindrical
ocean_stress::Vector{Float64}
atmosphere_stress::Vector{Float64}
+ thermal_energy::Float64
+
# Visualization parameters
color::Int
end
@@ -153,6 +155,8 @@ mutable struct GrainArrays
ocean_stress::Array{Float64, 2}
atmosphere_stress::Array{Float64, 2}
+ thermal_energy::Vector{Float64}
+
color::Vector{Int}
end
diff --git a/src/grain.jl b/src/grain.jl
@@ -25,7 +25,9 @@ export addGrainCylindrical!
rotating, enabled, verbose,
ocean_grid_pos, atmosphere_grid_pos,
n_contact, granular_stress, ocean_stress,
- atmosphere_stress])
+ atmosphere_stress,
+ thermal_energy,
+ color])
Creates and adds a cylindrical grain to a simulation. Most of the arguments
are optional, and come with default values. The only required arguments are
@@ -99,6 +101,7 @@ are optional, and come with default values. The only required arguments are
ocean drag [Pa].
* `atmosphere_stress::Vector{Float64} = [0., 0.]`: resultant stress on grain
from atmosphere drag [Pa].
+* `thermal_energy::Float64 = 0.0`: thermal energy of grain [J].
* `color::Int=0`: type number, usually used for associating a color to the grain
during visualization.
@@ -170,6 +173,7 @@ function addGrainCylindrical!(simulation::Simulation,
granular_stress::Vector{Float64} = [0., 0.],
ocean_stress::Vector{Float64} = [0., 0.],
atmosphere_stress::Vector{Float64} = [0., 0.],
+ thermal_energy::Float64 = 0.,
color::Int = 0)
# Check input values
@@ -267,6 +271,8 @@ function addGrainCylindrical!(simulation::Simulation,
ocean_stress,
atmosphere_stress,
+ thermal_energy,
+
color
)
@@ -435,6 +441,9 @@ function convertGrainDataToArrays(simulation::Simulation)
## atmosphere_stress
zeros(Float64, 3, length(simulation.grains)),
+ ## thermal_energy
+ Array{Float64}(length(simulation.grains)),
+
# Visualization parameters
## color
Array{Int}(length(simulation.grains)),
@@ -512,6 +521,8 @@ function convertGrainDataToArrays(simulation::Simulation)
ifarr.atmosphere_stress[1:2, i] =
simulation.grains[i].atmosphere_stress
+ ifarr.thermal_energy[i] = simulation.grains[i].thermal_energy
+
ifarr.color[i] = simulation.grains[i].color
end
@@ -578,6 +589,8 @@ function deleteGrainArrays!(ifarr::GrainArrays)
ifarr.ocean_stress = f2
ifarr.atmosphere_stress = f2
+ ifarr.thermal_energy = f1
+
ifarr.color = i1
gc()
nothing
@@ -643,7 +656,9 @@ function printGrainInfo(f::GrainCylindrical)
println(" granular_stress: $(f.granular_stress) Pa")
println(" ocean_stress: $(f.ocean_stress) Pa")
- println(" atmosphere_stress: $(f.atmosphere_stress) Pa")
+ println(" atmosphere_stress: $(f.atmosphere_stress) Pa\n")
+
+ println(" thermal_energy: $(f.thermal_energy) J\n")
println(" color: $(f.color)\n")
nothing
@@ -690,6 +705,26 @@ function totalGrainKineticRotationalEnergy(simulation::Simulation)
return E_sum
end
+export grainThermalEnergy
+"Returns the thermal energy of the grain, produced by Coulomb slip"
+function grainThermalEnergy(grain::GrainCylindrical)
+ return grain.thermal_energy
+end
+
+export totalGrainThermalEnergy
+"""
+ totalGrainKineticTranslationalEnergy(simulation)
+
+Returns the sum of thermal energy of all grains in a simulation
+"""
+function totalGrainThermalEnergy(simulation::Simulation)
+ E_sum = 0.
+ for grain in simulation.grains
+ E_sum += grainThermalEnergy(grain)
+ end
+ return E_sum
+end
+
export addBodyForce!
"""
setBodyForce!(grain, force)
@@ -791,6 +826,8 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.ocean_stress ≈ if2.ocean_stress
@test if1.atmosphere_stress ≈ if2.atmosphere_stress
+ @test if1.thermal_energy ≈ if2.thermal_energy
+
@test if1.color ≈ if2.color
nothing
end
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -256,9 +256,14 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
elseif k_t ≈ 0. && γ_t > 0.
force_t = abs(γ_t * vel_t)
+
+ # Coulomb slip
if force_t > μ_d_minimum*abs(force_n)
force_t = μ_d_minimum*abs(force_n)
- simulation.grains[i].contact_age[ic] = -simulation.time_step
+ simulation.grains[i].contact_age[ic] = 0.0
+ E_shear = abs(force_t)*vel_t*simulation.time_step
+ simulation.grains[i].thermal_energy += 0.5*E_shear
+ simulation.grains[j].thermal_energy += 0.5*E_shear
end
if vel_t > 0.
force_t = -force_t
@@ -268,10 +273,14 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
force_t = -k_t*δ_t - γ_t*vel_t
+ # Coulomb slip
if abs(force_t) > μ_d_minimum*abs(force_n)
force_t = μ_d_minimum*abs(force_n)*force_t/abs(force_t)
δ_t = (-force_t - γ_t*vel_t)/k_t
- simulation.grains[i].contact_age[ic] = -simulation.time_step
+ simulation.grains[i].contact_age[ic] = 0.0
+ E_shear = abs(force_t)*vel_t*simulation.time_step
+ simulation.grains[i].thermal_energy += 0.5*E_shear
+ simulation.grains[j].thermal_energy += 0.5*E_shear
end
else
diff --git a/src/io.jl b/src/io.jl
@@ -428,6 +428,9 @@ function writeGrainVTK(simulation::Simulation,
WriteVTK.vtk_point_data(vtkfile, ifarr.atmosphere_stress,
"Atmosphere stress [Pa]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.thermal_energy,
+ "Thermal energy [J]")
+
WriteVTK.vtk_point_data(vtkfile, ifarr.color,
"Color [-]")
diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
@@ -35,7 +35,11 @@ Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbos
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol
+E_thermal_final = Granular.totalGrainThermalEnergy(sim)
+println(E_kin_lin_init)
+println(E_kin_lin_final)
+println(E_thermal_final)
+@test E_kin_lin_init ≈ E_kin_lin_final+E_thermal_final atol=E_kin_lin_init*tol
@test E_kin_rot_init ≈ E_kin_rot_final
@@ -48,7 +52,8 @@ Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbos
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol
+E_thermal_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final+E_thermal_final atol=E_kin_lin_init*tol
@test E_kin_rot_init ≈ E_kin_rot_final
@@ -61,7 +66,8 @@ Granular.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verb
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol
+E_thermal_final = Granular.totalGrainThermalEnergy(sim)
+@test E_kin_lin_init ≈ E_kin_lin_final+E_thermal_final atol=E_kin_lin_init*tol
@test E_kin_rot_init ≈ E_kin_rot_final
info("# Ice floes free to move")
diff --git a/test/vtk.jl b/test/vtk.jl
@@ -29,7 +29,7 @@ end
grainpath = "test/test.grains.1.vtu"
grainchecksum =
-"8ce4e4045d150fc75cdb4126edcb32cb4d09f297065602a075cbe34dceef27c7 " *
+"f8aa4ab923b2466d7a6ffb835bfff4d0e0ff0a168ff267431a082bb3018caaa1 " *
grainpath * "\n"
graininteractionpath = "test/test.grain-interaction.1.vtp"
