ttemporal_integration.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
---
ttemporal_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