twall.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
twall.jl (9060B)
---
1 ## Manage dynamic walls in the model
2
3 export addWallLinearFrictionless!
4 """
5 function addWallLinear!(simulation, normal, pos[, bc, mass, thickness,
6 normal_stress, vel, acc, force,
7 contact_viscosity_normal, verbose])
8
9 Creates and adds a linear (flat) and frictionless dynamic wall to a grain to a
10 simulation. Most of the arguments are optional, and come with default values.
11 The only required arguments are
12 `simulation`, `normal`, `pos`, and `bc`.
13
14 # Arguments
15 * `simulation::Simulation`: the simulation object where the wall should be
16 added to.
17 * `normal::Vector{Float64}`: 3d vector denoting the normal to the wall [m]. The
18 wall will only interact in the opposite direction of this vector, so the
19 normal vector should point in the direction of the grains. If a 2d vector is
20 passed, the third (z) component is set to zero.
21 * `pos::Float64`: position along axis parallel to the normal vector [m].
22 * `bc::String="fixed"`: boundary condition, possible values are `"fixed"`
23 (default), `"normal stress"`, or `"velocity"`.
24 * `mass::Float64=NaN`: wall mass, which is used if wall boundary conditions
25 differs from `bc="fixed"`. If the parameter is left to its default value,
26 the wall mass is set to be equal the total mass of grains in the simulation.
27 Units: [kg]
28 * `thickness::Float64=NaN`: wall thickness, which is used for determining wall
29 surface area. If the parameter is left to its default value, the wall
30 thickness is set to be equal to the thickest grain in the simulation.
31 Units: [m].
32 * `normal_stress::Float64=0.`: the wall normal stress when `bc == "normal
33 stress"` [Pa].
34 * `vel::Float64=0.`: the wall velocity along the `normal` vector. If the
35 wall boundary condition is `bc = "velocity"` the wall will move according to
36 this constant value. If `bc = "normal stress"` the velocity will be a free
37 parameter. Units: [m/s]
38 * `force::Float64=0.`: sum of normal forces on the wall from interaction with
39 grains [N].
40 * `contact_viscosity_normal::Float64=0.`: viscosity to apply in parallel to
41 elasticity in interactions between wall and particles [N/(m/s)]. When this
42 term is larger than zero, the wall-grain interaction acts like a sink of
43 kinetic energy.
44 * `verbose::Bool=true`: show verbose information during function call.
45
46 # Examples
47 The most basic example adds a new fixed wall to the simulation `sim`, with a
48 wall-face normal of `[1., 0.]` (wall along *y* and normal to *x*), a position of
49 `1.5` meter:
50
51 ```julia
52 Granular.addWallLinearFrictionless!(sim, [1., 0., 0.], 1.5)
53 ```
54
55 The following example creates a wall with a velocity of 0.5 m/s towards *-y*:
56
57 ```julia
58 Granular.addWallLinearFrictionless!(sim, [0., 1., 0.], 1.5,
59 bc="velocity",
60 vel=-0.5)
61 ```
62
63 To create a wall parallel to the *y* axis pushing downwards with a constant
64 normal stress of 100 kPa, starting at a position of y = 3.5 m:
65
66 ```julia
67 Granular.addWallLinearFrictionless!(sim, [0., 1., 0.], 3.5,
68 bc="normal stress",
69 normal_stress=100e3)
70 ```
71 """
72 function addWallLinearFrictionless!(simulation::Simulation,
73 normal::Vector{Float64},
74 pos::Float64;
75 bc::String = "fixed",
76 mass::Float64 = -1.,
77 thickness::Float64 = -1.,
78 surface_area::Float64 = -1.,
79 normal_stress::Float64 = 0.,
80 vel::Float64 = 0.,
81 acc::Float64 = 0.,
82 force::Float64 = 0.,
83 contact_viscosity_normal::Float64 = 0.,
84 verbose::Bool=true)
85
86 # Check input values
87 if length(normal) != 3
88 normal = vecTo3d(normal)
89 end
90
91 if bc != "fixed" && bc != "velocity" && bc != "normal stress"
92 error("Wall BC must be 'fixed', 'velocity', or 'normal stress'.")
93 end
94
95 if !(normal ≈ [1., 0., 0.]) && !(normal ≈ [0., 1., 0.])
96 error("Currently only walls with normals orthogonal to the " *
97 "coordinate system are allowed, i.e. normals parallel to the " *
98 "x or y axes. Accepted values for `normal` " *
99 "are [1., 0., 0.] and [0., 1., 0.]. The passed normal was $normal")
100 end
101
102 # if not set, set wall mass to equal the mass of all grains.
103 if mass < 0.
104 if length(simulation.grains) < 1
105 error("If wall mass is not specified, walls should be " *
106 "added after grains have been added to the simulation.")
107 end
108 mass = 0.
109 for grain in simulation.grains
110 mass += grain.mass
111 end
112 if verbose
113 @info "Setting wall mass to total grain mass: $mass kg"
114 end
115 end
116
117 # if not set, set wall thickness to equal largest grain thickness
118 if thickness < 0.
119 if length(simulation.grains) < 1
120 error("If wall thickness is not specified, walls should " *
121 "be added after grains have been added to the simulation.")
122 end
123 thickness = -Inf
124 for grain in simulation.grains
125 if grain.thickness > thickness
126 thickness = grain.thickness
127 end
128 end
129 if verbose
130 @info "Setting wall thickness to max grain thickness: $thickness m"
131 end
132 end
133
134 # if not set, set wall surface area from the ocean grid
135 if surface_area < 0. && bc != "fixed"
136 if typeof(simulation.ocean.input_file) == Bool
137 error("simulation.ocean must be set beforehand")
138 end
139 surface_area = getWallSurfaceArea(simulation, normal, thickness)
140 end
141
142 # Create wall object
143 wall = WallLinearFrictionless(normal,
144 pos,
145 bc,
146 mass,
147 thickness,
148 surface_area,
149 normal_stress,
150 vel,
151 acc,
152 force,
153 contact_viscosity_normal)
154
155 # Add to simulation object
156 addWall!(simulation, wall, verbose)
157 nothing
158 end
159
160 export getWallSurfaceArea
161 """
162 getWallSurfaceArea(simulation, wall_index)
163
164 Returns the surface area of the wall given the grid size and its index.
165
166 # Arguments
167 * `simulation::Simulation`: the simulation object containing the wall.
168 * `wall_index::Integer=1`: the wall number in the simulation object.
169 """
170 function getWallSurfaceArea(sim::Simulation, wall_index::Integer)
171
172 if sim.walls[wall_index].normal ≈ [1., 0., 0.]
173 return (sim.ocean.yq[end,end] - sim.ocean.yq[1,1]) *
174 sim.walls[wall_index].thickness
175 elseif sim.walls[wall_index].normal ≈ [0., 1., 0.]
176 return (sim.ocean.xq[end,end] - sim.ocean.xq[1,1]) *
177 sim.walls[wall_index].thickness
178 else
179 error("Wall normal not understood")
180 end
181 nothing
182 end
183 function getWallSurfaceArea(sim::Simulation, normal::Vector{Float64},
184 thickness::Float64)
185
186 if length(normal) == 3 && normal ≈ [1., 0., 0.]
187 return (sim.ocean.yq[end,end] - sim.ocean.yq[1,1]) * thickness
188 elseif length(normal) == 3 && normal ≈ [0., 1., 0.]
189 return (sim.ocean.xq[end,end] - sim.ocean.xq[1,1]) * thickness
190 else
191 error("Wall normal not understood")
192 end
193 nothing
194 end
195
196 export getWallNormalStress
197 """
198 getWallNormalStress(simulation[, wall_index, stress_type])
199
200 Returns the current "effective" or "defined" normal stress on the wall with
201 index `wall_index` inside the `simulation` object. The returned value is given
202 in Pascal.
203
204 # Arguments
205 * `simulation::Simulation`: the simulation object containing the wall.
206 * `wall_index::Integer=1`: the wall number in the simulation object.
207 * `stress_type::String="effective"`: the normal-stress type to return. The
208 defined value corresponds to the normal stress that the wall is asked to
209 uphold. The effective value is the actual current normal stress. Usually,
210 the magnitude of the effective normal stress fluctuates around the defined
211 normal stress.
212 """
213 function getWallNormalStress(sim::Simulation;
214 wall_index::Integer=1,
215 stress_type::String="effective")
216 if stress_type == "defined"
217 return sim.walls[wall_index].normal_stress
218
219 elseif stress_type == "effective"
220 return sim.walls[wall_index].force / getWallSurfaceArea(sim, wall_index)
221 else
222 error("stress_type not understood, " *
223 "should be 'effective' or 'defined' but is '$stress_type'.")
224 end
225 end