temporal_integration.jl (8634B)
1 2 export updateGrainKinematics! 3 """ 4 updateGrainKinematics!(simulation::Simulation[, 5 method::String = "Three-term Taylor"]) 6 7 Update the grain kinematic parameters using a temporal integration scheme, 8 the current force and torque balance, and gravitational acceleration. If the 9 simulation contains a grid with periodic boundaries, affected grain positions 10 are adjusted accordingly. 11 12 # Arguments 13 * `simulation::Simulation`: update the grain positions in this object 14 according to temporal integration of length `simulation.time_step`. 15 * `method::String = "Three-term Taylor"`: the integration method to use. 16 Available methods are "Two-term Taylor" and "Three-term Taylor". The 17 three-term Taylor expansion is slightly more computationally expensive than 18 the two-term Taylor expansion, but offers an order-of-magnitude increase in 19 precision of grain positions. The two-term expansion can obtain similar 20 precision if the time step is 1/10 the length. 21 """ 22 function updateGrainKinematics!(simulation::Simulation; 23 method::String = "Three-term Taylor") 24 25 if method == "Two-term Taylor" 26 for grain in simulation.grains 27 if !grain.enabled 28 continue 29 end 30 updateGrainKinematicsTwoTermTaylor!(grain, simulation) 31 end 32 elseif method == "Three-term Taylor" 33 for grain in simulation.grains 34 if !grain.enabled 35 continue 36 end 37 updateGrainKinematicsThreeTermTaylor!(grain, simulation) 38 end 39 else 40 error("Unknown integration method '$method'") 41 end 42 moveGrainsAcrossPeriodicBoundaries!(simulation) 43 reflectGrainsFromImpermeableBoundaries!(simulation) 44 nothing 45 end 46 47 export updateGrainKinematicsTwoTermTaylor! 48 """ 49 Use a two-term Taylor expansion for integrating the kinematic degrees of freedom 50 for a `grain`. 51 """ 52 function updateGrainKinematicsTwoTermTaylor!(grain::GrainCylindrical, 53 simulation::Simulation) 54 grain.lin_acc = grain.force/grain.mass 55 grain.ang_acc = grain.torque/grain.moment_of_inertia 56 57 if grain.fixed 58 if !grain.allow_x_acc 59 grain.lin_acc[1] = 0. 60 end 61 if !grain.allow_y_acc 62 grain.lin_acc[2] = 0. 63 end 64 if !grain.allow_z_acc 65 grain.lin_acc[3] = 0. 66 end 67 grain.ang_acc = zeros(3) 68 elseif !grain.rotating 69 grain.ang_acc = zeros(3) 70 end 71 72 grain.lin_pos += 73 grain.lin_vel * simulation.time_step + 74 0.5*grain.lin_acc * simulation.time_step^2.0 75 76 grain.lin_disp += 77 grain.lin_vel * simulation.time_step + 78 0.5*grain.lin_acc * simulation.time_step^2.0 79 80 grain.ang_pos += 81 grain.ang_vel * simulation.time_step + 82 0.5*grain.ang_acc * simulation.time_step^2.0 83 84 grain.lin_vel += grain.lin_acc * simulation.time_step 85 grain.ang_vel += grain.ang_acc * simulation.time_step 86 nothing 87 end 88 89 export updateGrainKinematicsThreeTermTaylor! 90 """ 91 Use a three-term Taylor expansion for integrating the kinematic degrees of 92 freedom for a `grain`. 93 """ 94 function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical, 95 simulation::Simulation) 96 97 if simulation.time_iteration == 0 98 lin_acc_0 = zeros(3) 99 ang_acc_0 = zeros(3) 100 else 101 lin_acc_0 = grain.lin_acc 102 ang_acc_0 = grain.ang_acc 103 end 104 105 grain.lin_acc = grain.force/grain.mass 106 grain.ang_acc = grain.torque/grain.moment_of_inertia 107 108 if grain.fixed 109 if !grain.allow_x_acc 110 grain.lin_acc[1] = 0. 111 end 112 if !grain.allow_y_acc 113 grain.lin_acc[2] = 0. 114 end 115 if !grain.allow_z_acc 116 grain.lin_acc[3] = 0. 117 end 118 grain.ang_acc = zeros(3) 119 elseif !grain.rotating 120 grain.ang_acc = zeros(3) 121 end 122 123 # Temporal gradient in acceleration using backwards differences 124 d_lin_acc_dt = (grain.lin_acc - lin_acc_0)/simulation.time_step 125 d_ang_acc_dt = (grain.ang_acc - ang_acc_0)/simulation.time_step 126 127 grain.lin_pos += 128 grain.lin_vel * simulation.time_step + 129 0.5 * grain.lin_acc * simulation.time_step^2. + 130 1. / 6. * d_lin_acc_dt * simulation.time_step^3. 131 132 grain.lin_disp += 133 grain.lin_vel * simulation.time_step + 134 0.5 * grain.lin_acc * simulation.time_step^2. + 135 1. / 6. * d_lin_acc_dt * simulation.time_step^3. 136 137 grain.ang_pos += 138 grain.ang_vel * simulation.time_step + 139 0.5 * grain.ang_acc * simulation.time_step^2. + 140 1. / 6. * d_ang_acc_dt * simulation.time_step^3. 141 142 grain.lin_vel += grain.lin_acc * simulation.time_step + 143 0.5 * d_lin_acc_dt * simulation.time_step^2. 144 grain.ang_vel += grain.ang_acc * simulation.time_step + 145 0.5 * d_ang_acc_dt * simulation.time_step^2. 146 nothing 147 end 148 149 export updateWallKinematics! 150 """ 151 updateWallKinematics!(simulation::Simulation[, 152 method::String = "Three-term Taylor"]) 153 154 Update the wall kinematic parameters using a temporal integration scheme, 155 the current force and torque balance, and gravitational acceleration. If the 156 simulation contains a grid with periodic boundaries, affected wall positions 157 are adjusted accordingly. 158 159 # Arguments 160 * `simulation::Simulation`: update the wall positions in this object 161 according to temporal integration of length `simulation.time_step`. 162 * `method::String = "Three-term Taylor"`: the integration method to use. 163 Available methods are "Two-term Taylor" and "Three-term Taylor". The 164 three-term Taylor expansion is slightly more computationally expensive than 165 the two-term Taylor expansion, but offers an order-of-magnitude increase in 166 precision of wall positions. The two-term expansion can obtain similar 167 precision if the time step is 1/10 the length. 168 """ 169 function updateWallKinematics!(simulation::Simulation; 170 method::String = "Three-term Taylor") 171 172 if method == "Two-term Taylor" 173 for wall in simulation.walls 174 if wall.bc == "fixed" 175 continue 176 end 177 updateWallKinematicsTwoTermTaylor!(wall, simulation) 178 end 179 elseif method == "Three-term Taylor" 180 for wall in simulation.walls 181 if wall.bc == "fixed" 182 continue 183 end 184 updateWallKinematicsThreeTermTaylor!(wall, simulation) 185 end 186 else 187 error("Unknown integration method '$method'") 188 end 189 nothing 190 end 191 192 export updateWallKinematicsTwoTermTaylor! 193 """ 194 updateWallKinematicsTwoTermTaylor!(wall, simulation) 195 196 Use a two-term Taylor expansion for integrating the kinematic degrees of freedom 197 for a `wall`. 198 """ 199 function updateWallKinematicsTwoTermTaylor!(wall::WallLinearFrictionless, 200 simulation::Simulation) 201 if wall.bc == "fixed" 202 return nothing 203 elseif wall.bc == "velocity" 204 wall.acc = 0.0 205 elseif wall.bc == "normal stress" 206 wall.acc = (wall.force + wall.normal_stress*wall.surface_area)/wall.mass 207 else 208 error("wall boundary condition was not understood ($(wall.bc))") 209 end 210 211 wall.pos += 212 wall.vel * simulation.time_step + 213 0.5*wall.acc * simulation.time_step^2.0 214 215 wall.vel += wall.acc * simulation.time_step 216 nothing 217 end 218 219 220 export updateWallKinematicsThreeTermTaylor! 221 """ 222 updateWallKinematicsThreeTermTaylor!(wall, simulation) 223 224 Use a two-term Taylor expansion for integrating the kinematic degrees of freedom 225 for a `wall`. 226 """ 227 function updateWallKinematicsThreeTermTaylor!(wall::WallLinearFrictionless, 228 simulation::Simulation) 229 if simulation.time_iteration == 0 230 acc_0 = 0. 231 else 232 acc_0 = wall.acc 233 end 234 235 if wall.bc == "fixed" 236 return nothing 237 elseif wall.bc == "velocity" 238 wall.acc = 0.0 239 elseif wall.bc == "normal stress" 240 wall.acc = (wall.force + wall.normal_stress*wall.surface_area)/wall.mass 241 else 242 error("wall boundary condition was not understood ($(wall.bc))") 243 end 244 245 # Temporal gradient in acceleration using backwards differences 246 d_acc_dt = (wall.acc - acc_0)/simulation.time_step 247 248 wall.pos += 249 wall.vel * simulation.time_step + 250 0.5*wall.acc * simulation.time_step^2.0 + 251 1.0/6.0 * d_acc_dt * simulation.time_step^3.0 252 253 wall.vel += wall.acc * simulation.time_step + 254 0.5 * d_acc_dt * simulation.time_step^2.0 255 256 nothing 257 end