wall.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