tcfd.rst - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tcfd.rst (21635B)
---
1 Fluid simulation and particle-fluid interaction
2 ===============================================
3 A new and experimental addition to *sphere* is the ability to simulate a mixture
4 of particles and a Newtonian fluid. The fluid is simulated using an Eulerian
5 continuum approach, using a custom CUDA solver for GPU computation. This
6 approach allows for fast simulations due to the limited need for GPU-CPU
7 communications, as well as a flexible code base.
8
9 The following sections will describe the theoretical background, as well as the
10 solution procedure and the numerical implementation.
11
12 Derivation of the Navier Stokes equations with porosity
13 -------------------------------------------------------
14 Following the outline presented by `Limache and Idelsohn (2006)`_, the
15 continuity equation for an incompressible fluid material is given by:
16
17 .. math::
18 \nabla \cdot \boldsymbol{v} = 0
19
20 and the momentum equation:
21
22 .. math::
23 \rho \frac{\partial \boldsymbol{v}}{\partial t}
24 + \rho (\boldsymbol{v} \cdot \nabla \boldsymbol{v})
25 = \nabla \cdot \boldsymbol{\sigma}
26 - \boldsymbol{f}^i
27 + \rho \boldsymbol{g}
28
29 Here, :math:`\boldsymbol{v}` is the fluid velocity, :math:`\rho` is the
30 fluid density, :math:`\boldsymbol{\sigma}` is the `Cauchy stress tensor`_,
31 :math:`\boldsymbol{f}^i` is the particle-fluid interaction vector and
32 :math:`\boldsymbol{g}` is the gravitational acceleration. For incompressible
33 Newtonian fluids, the Cauchy stress is given by:
34
35 .. math::
36 \boldsymbol{\sigma} = -p \boldsymbol{I} + \boldsymbol{\tau}
37
38 :math:`p` is the fluid pressure, :math:`\boldsymbol{I}` is the identity
39 tensor, and :math:`\boldsymbol{\tau}` is the deviatoric stress tensor, given
40 by:
41
42 .. math::
43 \boldsymbol{\tau} =
44 \mu_f \nabla \boldsymbol{v}
45 + \mu_f (\nabla \boldsymbol{v})^T
46
47 By using the following vector identities:
48
49 .. math::
50 \nabla \cdot (p \boldsymbol{I}) = \nabla p
51
52 \nabla \cdot (\nabla \boldsymbol{v}) = \nabla^2 \boldsymbol{v}
53
54 \nabla \cdot (\nabla \boldsymbol{v})^T
55 = \nabla (\nabla \cdot \boldsymbol{v})
56
57 the deviatoric component of the Cauchy stress tensor simplifies to the
58 following, assuming that spatial variations in the viscosity can be neglected:
59
60 .. math::
61 = -\nabla p
62 + \mu_f \nabla^2 \boldsymbol{v}
63
64 Since we are dealing with fluid flow in a porous medium, additional terms are
65 introduced to the equations for conservation of mass and momentum. In the
66 following, the equations are derived for the first spatial component. The
67 solution for the other components is trivial.
68
69 The porosity value (in the saturated porous medium the volumetric fraction of
70 the fluid phase) denoted :math:`\phi` is incorporated in the continuity and
71 momentum equations. The continuity equation becomes:
72
73 .. math::
74 \frac{\partial \phi}{\partial t}
75 + \nabla \cdot (\phi \boldsymbol{v}) = 0
76
77 For the :math:`x` component, the Lagrangian formulation of the momentum equation
78 with a body force :math:`\boldsymbol{f}` becomes:
79
80 .. math::
81 \frac{D (\phi v_x)}{D t}
82 = \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\sigma}) \right]_x
83 - \frac{1}{\rho} f^i_x
84 + \phi g
85
86 In the Eulerian formulation, an advection term is added, and the Cauchy stress
87 tensor is represented as isotropic and deviatoric components individually:
88
89 .. math::
90 \frac{\partial (\phi v_x)}{\partial t}
91 + \boldsymbol{v} \cdot \nabla (\phi v_x)
92 = \frac{1}{\rho} \left[ \nabla \cdot (-\phi p \boldsymbol{I})
93 + \phi \boldsymbol{\tau}) \right]_x
94 - \frac{1}{\rho} f^i_x
95 + \phi g_x
96
97 Using vector identities to rewrite the advection term, and expanding the fluid
98 stress tensor term:
99
100 .. math::
101 \frac{\partial (\phi v_x)}{\partial t}
102 + \nabla \cdot (\phi v_x \boldsymbol{v})
103 - \phi v_x (\nabla \cdot \boldsymbol{v})
104 = \frac{1}{\rho} \left[ -\nabla \phi p \right]_x
105 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
106 - \frac{1}{\rho} f^i_x
107 + \phi g_x
108
109 Spatial variations in the porosity are neglected,
110
111 .. math::
112 \nabla \phi := 0
113
114 and the pressure is attributed to the fluid phase alone (model B in Zhu et al.
115 2007 and Zhou et al. 2010). The divergence of fluid velocities is defined to be
116 zero:
117
118 .. math::
119 \nabla \cdot \boldsymbol{v} := 0
120
121 With these assumptions, the momentum equation simplifies to:
122
123 .. math::
124 \frac{\partial (\phi v_x)}{\partial t}
125 + \nabla \cdot (\phi v_x \boldsymbol{v})
126 = -\frac{1}{\rho} \frac{\partial p}{\partial x}
127 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
128 - \frac{1}{\rho} f^i_x
129 + \phi g_x
130
131 The remaining part of the advection term is for the :math:`x` component
132 found as:
133
134 .. math::
135 \nabla \cdot (\phi v_x \boldsymbol{v}) =
136 \left[
137 \frac{\partial}{\partial x},
138 \frac{\partial}{\partial y},
139 \frac{\partial}{\partial z}
140 \right]
141 \left[
142 \begin{array}{c}
143 \phi v_x v_x\\
144 \phi v_x v_y\\
145 \phi v_x v_z\\
146 \end{array}
147 \right]
148 =
149 \frac{\partial (\phi v_x v_x)}{\partial x} +
150 \frac{\partial (\phi v_x v_y)}{\partial y} +
151 \frac{\partial (\phi v_x v_z)}{\partial z}
152
153 The deviatoric stress tensor is in this case symmetrical, i.e. :math:`\tau_{ij}
154 = \tau_{ji}`, and is found by:
155
156 .. math::
157 \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
158 = \frac{1}{\rho}
159 \left[
160 \left[
161 \frac{\partial}{\partial x},
162 \frac{\partial}{\partial y},
163 \frac{\partial}{\partial z}
164 \right]
165 \phi
166 \left[
167 \begin{matrix}
168 \tau_{xx} & \tau_{xy} & \tau_{xz}\\
169 \tau_{yx} & \tau_{yy} & \tau_{yz}\\
170 \tau_{zx} & \tau_{zy} & \tau_{zz}\\
171 \end{matrix}
172 \right]
173 \right]_x
174
175 = \frac{1}{\rho}
176 \left[
177 \begin{array}{c}
178 \frac{\partial (\phi \tau_{xx})}{\partial x}
179 + \frac{\partial (\phi \tau_{xy})}{\partial y}
180 + \frac{\partial (\phi \tau_{xz})}{\partial z}\\
181 \frac{\partial (\phi \tau_{yx})}{\partial x}
182 + \frac{\partial (\phi \tau_{yy})}{\partial y}
183 + \frac{\partial (\phi \tau_{yz})}{\partial z}\\
184 \frac{\partial (\phi \tau_{zx})}{\partial x}
185 + \frac{\partial (\phi \tau_{zy})}{\partial y}
186 + \frac{\partial (\phi \tau_{zz})}{\partial z}\\
187 \end{array}
188 \right]_x
189 = \frac{1}{\rho}
190 \left(
191 \frac{\partial (\phi \tau_{xx})}{\partial x}
192 + \frac{\partial (\phi \tau_{xy})}{\partial y}
193 + \frac{\partial (\phi \tau_{xz})}{\partial z}
194 \right)
195
196 In a linear viscous fluid, the stress and strain rate
197 (:math:`\dot{\boldsymbol{\epsilon}}`) is linearly dependent, scaled by the
198 viscosity parameter :math:`\mu_f`:
199
200 .. math::
201 \tau_{ij} = 2 \mu_f \dot{\epsilon}_{ij}
202 = \mu_f \left(
203 \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}
204 \right)
205
206 With this relationship, the deviatoric stress tensor components can be
207 calculated as:
208
209 .. math::
210 \tau_{xx} = 2 \mu_f \frac{\partial v_x}{\partial x} \qquad
211 \tau_{yy} = 2 \mu_f \frac{\partial v_y}{\partial y} \qquad
212 \tau_{zz} = 2 \mu_f \frac{\partial v_z}{\partial z}
213
214 \tau_{xy} = \mu_f \left(
215 \frac{\partial v_x}{\partial y} + \frac{\partial v_y}{\partial x} \right)
216
217 \tau_{xz} = \mu_f \left(
218 \frac{\partial v_x}{\partial z} + \frac{\partial v_z}{\partial x} \right)
219
220 \tau_{yz} = \mu_f \left(
221 \frac{\partial v_y}{\partial z} + \frac{\partial v_z}{\partial y} \right)
222
223 where :math:`\mu_f` is the dynamic viscosity. The above formulation of the
224 fluid rheology assumes identical bulk and shear viscosities. The derivation of
225 the equations for the other spatial components is trivial.
226
227 Porosity estimation
228 -------------------
229 The solid volume in each fluid cell is determined by the ratio of the
230 a cell-centered spherical cell volume (:math:`V_c`) and the sum of intersecting
231 particle volumes (:math:`V_s`). The spherical cell volume has a center at
232 :math:`\boldsymbol{x}_i`, and a radius of :math:`R_i`, which is equal to half
233 the fluid cell width. The nearby particles are characterized by position
234 :math:`\boldsymbol{x}_j` and radius :math:`r_j`. The center distance is defined
235 as:
236
237 .. math::
238 d_{ij} = ||\boldsymbol{x}_i - \boldsymbol{x}_j||
239
240 The common volume of the two intersecting spheres is zero if the volumes aren't
241 intersecting, lens shaped if they are intersecting, and spherical if the
242 particle is fully contained by the spherical cell volume:
243
244 .. math::
245 V^s_{i} = \sum_j
246 \begin{cases}
247 0 & \textit{if } R_i + r_j \leq d_{ij} \\
248 \frac{1}{12d_{ij}} \left[ \pi (R_i + r_j - d_{ij})^2
249 (d_{ij}^2 + 2d_{ij}r_j - 3r_j^2 + 2d_{ij} R_i + 6r_j R_i - 3R_i^2)
250 \right] & \textit{if } R_i - r_j < d_{ij} < R_i + r_j \\
251 \frac{4}{3} \pi r^3_j & \textit{if } d_{ij} \leq R_i - r_j
252 \end{cases}
253
254 Using this method, the cell porosity values are continuous through time as
255 particles enter and exit the cell volume. The rate of porosity change
256 (:math:`d\phi/dt`) is estimated by the backwards Euler method
257 by considering the previous and current porosity.
258
259 Particle-fluid interaction
260 --------------------------
261 The momentum exchange of the granular and fluid phases follows the procedure
262 outlined by Gidaspow 1992 and Shamy and Zhegal 2005. The fluid and particle
263 interaction is based on the concept of drag, where the magnitude is based on
264 semi-empirical relationships. The drag force scales linearly with the relative
265 difference in velocity between the fluid and particle phase. On the base of
266 Newton's third law, the resulting drag force is applied with opposite signs to
267 the particle and fluid.
268
269 For fluid cells with porosities (:math:`\phi`) less or equal to 0.8, the drag
270 force is based on the Ergun (1952) equation:
271
272 .. math::
273 \bar{\boldsymbol{f}}_d = \left(
274 150 \frac{\mu_f (1-\phi)^2}{\phi\bar{d}^2}
275 + 1.75 \frac{(1-\phi)\rho_f
276 ||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||}{\bar{d}}
277 \right)
278 (\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p)
279
280 here, :math:`\bar{d}` denotes the average particle diameter in the cell,
281 :math:`\boldsymbol{v}_f` is the fluid flow velocity, and
282 :math:`\bar{\boldsymbol{v}}_p` is the average particle velocity in the cell. All
283 particles in contact with the previously mentioned cell-centered sphere for
284 porosity estimation contribute to the average particle velocity and diameter in
285 the fluid cell.
286
287 If the porosity is greater than 0.8, the cell-averaged drag force
288 (:math:`\bar{\boldsymbol{f}}_d` is found from the Wen and Yu (1966) equation,
289 which considers the fluid flow situation:
290
291 .. math::
292 \bar{\boldsymbol{f}}_d = \left(
293 \frac{3}{4}
294 \frac{C_d (1-\phi) \phi^{-2.65} \mu_f \rho_f
295 ||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||}{\bar{d}}
296 \right)
297 (\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p)
298
299 The drag coefficient :math:`C_d` is evaluated depending on the magnitude of the
300 Reynolds number :math:`Re`:
301
302 .. math::
303 C_d =
304 \begin{cases}
305 \frac{24}{Re} (1+0.15 (Re)^{0.687} & \textit{if } Re < 1,000 \\
306 0.44 & \textit{if } Re \geq 1,000
307 \end{cases}
308
309 where the Reynold's number is found by:
310
311 .. math::
312 Re = \frac{\phi\rho_f\bar{d}}{\mu_f}
313 ||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||
314
315 The interaction force is applied to the fluid with negative sign as a
316 contribution to the body force :math:`\boldsymbol{f}`. The fluid interaction
317 force applied particles in the fluid cell is:
318
319 .. math::
320 \boldsymbol{f}_i = \frac{\bar{\boldsymbol{f}}_d V_p}{1-\phi}
321
322 where :math:`V_p` denotes the particle volume. Optionally, the above
323 interaction force could be expanded to include the force induced by the fluid
324 pressure gradient:
325
326 .. math::
327 \boldsymbol{f}_i = \left(
328 -\nabla p +
329 \frac{\bar{\boldsymbol{f}}_d}{1-\phi}
330 \right) V_p
331
332
333 Fluid dynamics solution procedure by operator splitting
334 -------------------------------------------------------
335 The partial differential terms in the previously described equations are found
336 using finite central differences. Modifying the operator splitting methodology
337 presented by Langtangen et al. (2002), the predicted velocity
338 :math:`\boldsymbol{v}^*` after a finite time step
339 :math:`\Delta t` is found by explicit integration of the momentum equation.
340
341 .. math::
342 \frac{\Delta (\phi v_x)}{\Delta t}
343 + \nabla \cdot (\phi v_x \boldsymbol{v})
344 = - \frac{1}{\rho} \frac{\Delta p}{\Delta x}
345 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
346 - \frac{1}{\rho} f^i_x
347 + \phi g_x
348
349 \Downarrow
350
351 \phi \frac{\Delta v_x}{\Delta t}
352 + v_x \frac{\Delta \phi}{\Delta t}
353 + \nabla \cdot (\phi v_x \boldsymbol{v})
354 = - \frac{1}{\rho} \frac{\Delta p}{\Delta x}
355 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
356 - \frac{1}{\rho} f^i_x
357 + \phi g_x
358
359 We want to isolate :math:`\Delta v_x` in the above equation in order to project
360 the new velocity.
361
362 .. math::
363 \phi \frac{\Delta v_x}{\Delta t}
364 = - \frac{1}{\rho} \frac{\Delta p}{\Delta x}
365 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
366 - \frac{1}{\rho} f^i_x
367 + \phi g_x
368 - v_x \frac{\Delta \phi}{\Delta t}
369 - \nabla \cdot (\phi v_x \boldsymbol{v})
370
371 \Delta v_x
372 = - \frac{1}{\rho} \frac{\Delta p}{\Delta x} \frac{\Delta t}{\phi}
373 + \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
374 \frac{\Delta t}{\phi}
375 - \frac{\Delta t}{\rho\phi} f^i_x
376 + \Delta t g_x
377 - v_x \frac{\Delta \phi}{\phi}
378 - \nabla \cdot (\phi v_x \boldsymbol{v}) \frac{\Delta t}{\phi}
379
380 The term :math:`\beta` is introduced as an adjustable, dimensionless parameter
381 in the range :math:`[0;1]`, and determines the importance of the old pressure
382 values in the solution procedure (Langtangen et al. 2002). A value of 0
383 corresponds to `Chorin's projection method`_ originally described
384 in `Chorin (1968)`_.
385
386 .. math::
387 v_x^* = v_x^t + \Delta v_x
388
389 v_x^* = v_x^t
390 - \frac{\beta}{\rho} \frac{\Delta p^t}{\Delta x} \frac{\Delta t}{\phi^t}
391 + \frac{1}{\rho} \left[ \nabla \cdot (\phi^t \boldsymbol{\tau}^t) \right]_x
392 \frac{\Delta t}{\phi}
393 - \frac{\Delta t}{\rho\phi} f^i_x
394 + \Delta t g_x
395 - v^t_x \frac{\Delta \phi}{\phi^t}
396 - \nabla \cdot (\phi^t v_x^t \boldsymbol{v}^t) \frac{\Delta t}{\phi^t}
397
398 Here, :math:`\Delta x` denotes the cell spacing. The velocity found
399 (:math:`v_x^*`) is only a prediction of the fluid velocity at time
400 :math:`t+\Delta t`, since the estimate isn't constrained by the continuity
401 equation:
402
403 .. math::
404 \frac{\Delta \phi^t}{\Delta t} + \nabla \cdot (\phi^t
405 \boldsymbol{v}^{t+\Delta t}) = 0
406
407 The divergence of a scalar and vector can be `split`_:
408
409 .. math::
410 \phi^t \nabla \cdot \boldsymbol{v}^{t+\Delta t} +
411 \boldsymbol{v}^{t+\Delta t} \cdot \nabla \phi^t
412 + \frac{\Delta \phi^t}{\Delta t} = 0
413
414 The predicted velocity is corrected using the new pressure (Langtangen et al.
415 2002):
416
417 .. math::
418 \boldsymbol{v}^{t+\Delta t} = \boldsymbol{v}^*
419 %- \frac{\Delta t}{\rho} \nabla \epsilon
420 - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon
421 \quad \text{where} \quad
422 \epsilon = p^{t+\Delta t} - \beta p^t
423
424 The above formulation of the future velocity is put into the continuity
425 equation:
426
427 .. math::
428 \Rightarrow
429 \phi^t \nabla \cdot
430 \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right)
431 +
432 \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right)
433 \cdot \nabla \phi^t + \frac{\Delta \phi^t}{\Delta t} = 0
434
435 .. math::
436 \Rightarrow
437 \phi^t \nabla \cdot
438 \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \phi^t \nabla^2 \epsilon
439 + \nabla \phi^t \cdot \boldsymbol{v}^*
440 - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t}
441 + \frac{\Delta \phi^t}{\Delta t} = 0
442
443 .. math::
444 \Rightarrow
445 \frac{\Delta t}{\rho} \nabla^2 \epsilon
446 = \phi^t \nabla \cdot \boldsymbol{v}^*
447 + \nabla \phi^t \cdot \boldsymbol{v}^*
448 - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t}
449 + \frac{\Delta \phi^t}{\Delta t}
450
451 The pressure difference in time becomes a `Poisson equation`_ with added terms:
452
453 .. math::
454 \Rightarrow
455 \nabla^2 \epsilon
456 = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t}
457 + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t}
458 - \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t}
459 + \frac{\Delta \phi^t \rho}{\Delta t^2}
460
461 The right hand side of the above equation is termed the *forcing function*
462 :math:`f`, which is decomposed into two terms, :math:`f_1` and :math:`f_2`:
463
464 .. math::
465 f_1
466 = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t}
467 + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t}
468 + \frac{\Delta \phi^t \rho}{\Delta t^2}
469
470 f_2 =
471 \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t}
472
473
474 During the `Jacobi iterative solution procedure`_ :math:`f_1` remains constant,
475 while :math:`f_2` changes value. For this reason, :math:`f_1` is found only
476 during the first iteration, while :math:`f_2` is updated every time. The value
477 of the forcing function is found as:
478
479 .. math::
480 f = f_1 - f_2
481
482 Using second-order finite difference approximations of the Laplace operator
483 second-order partial derivatives, the differential equations become a system of
484 equations that is solved using `iteratively`_ using Jacobi updates. The total
485 number of unknowns is :math:`(n_x - 1)(n_y - 1)(n_z - 1)`.
486
487 The discrete Laplacian (approximation of the Laplace operator) can be obtained
488 by a finite-difference seven-point stencil in a three-dimensional, cubic
489 grid with cell spacing :math:`\Delta x, \Delta y, \Delta z`, considering the six
490 face neighbors:
491
492 .. math::
493 \nabla^2 \epsilon_{i_x,i_y,i_z} \approx
494 \frac{\epsilon_{i_x-1,i_y,i_z} - 2 \epsilon_{i_x,i_y,i_z}
495 + \epsilon_{i_x+1,i_y,i_z}}{\Delta x^2}
496 + \frac{\epsilon_{i_x,i_y-1,i_z} - 2 \epsilon_{i_x,i_y,i_z}
497 + \epsilon_{i_x,i_y+1,i_z}}{\Delta y^2}
498
499 + \frac{\epsilon_{i_x,i_y,i_z-1} - 2 \epsilon_{i_x,i_y,i_z}
500 + \epsilon_{i_x,i_y,i_z+1}}{\Delta z^2}
501 \approx f_{i_x,i_y,i_z}
502
503 Within a Jacobi iteration, the value of the unknowns (:math:`\epsilon^n`) is
504 used to find an updated solution estimate (:math:`\epsilon^{n+1}`).
505 The solution for the updated value takes the form:
506
507 .. math::
508 \epsilon^{n+1}_{i_x,i_y,i_z}
509 = \frac{-\Delta x^2 \Delta y^2 \Delta z^2 f_{i_x,i_y,i_z}
510 + \Delta y^2 \Delta z^2 (\epsilon^n_{i_x-1,i_y,i_z} +
511 \epsilon^n_{i_x+1,i_y,i_z})
512 + \Delta x^2 \Delta z^2 (\epsilon^n_{i_x,i_y-1,i_z} +
513 \epsilon^n_{i_x,i_y+1,i_z})
514 + \Delta x^2 \Delta y^2 (\epsilon^n_{i_x,i_y,i_z-1} +
515 \epsilon^n_{i_x,i_y,i_z+1})}
516 {2 (\Delta x^2 \Delta y^2
517 + \Delta x^2 \Delta z^2
518 + \Delta y^2 \Delta z^2) }
519
520 The difference between the current and updated value is termed the *normalized
521 residual*:
522
523 .. math::
524 r_{i_x,i_y,i_z} = \frac{(\epsilon^{n+1}_{i_x,i_y,i_z}
525 - \epsilon^n_{i_x,i_y,i_z})^2}{(\epsilon^{n+1}_{i_x,i_y,i_z})^2}
526
527 Note that the :math:`\epsilon` values cannot be 0 due to the above normalization
528 of the residual.
529
530 The updated values are at the end of the iteration stored as the current values,
531 and the maximal value of the normalized residual is found. If this value is
532 larger than a tolerance criteria, the procedure is repeated. The iterative
533 procedure is ended if the number of iterations exceeds a defined limit.
534
535 After the values of :math:`\epsilon` are found, they are used to find the new
536 pressures and velocities:
537
538 .. math::
539 \bar{p}^{t+\Delta t} = \beta \bar{p}^t + \epsilon
540
541 .. math::
542 \bar{\boldsymbol{v}}^{t+\Delta t} =
543 \bar{\boldsymbol{v}}^* - \frac{\Delta t}{\rho\phi} \nabla \epsilon
544
545
546 Boundary conditions
547 -------------------
548 The lateral boundaries are periodic. This cannot be changed in the current
549 version of ``sphere``. This means that the fluid properties at the paired,
550 parallel lateral (:math:`x` and :math:`y`) boundaries are identical. A flow
551 leaving through one side reappears on the opposite side.
552
553 The top and bottom boundary conditions of the fluid grid can be either:
554 prescribed pressure (Dirichlet), or prescribed velocity (Neumann). The
555 (horizontal) velocities parallel to the boundaries are free to attain other
556 values (free slip). The Dirichlet boundary condition is enforced by keeping the
557 value of :math:`\epsilon` constant at the boundaries, e.g.:
558
559 .. math::
560 \epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z}
561 =
562 \epsilon^{n}_{i_x,i_y,i_z = 1 \vee n_z}
563
564 The Neumann boundary condition of no flow across the boundary is enforced by
565 setting the gradient of :math:`\epsilon` perpendicular to the boundary to zero,
566 e.g.:
567
568 .. math::
569 \nabla_z \epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z} = 0
570
571
572 Numerical implementation
573 ------------------------
574 Ghost nodes
575
576 ---
577
578
579
580
581 .. _Limache and Idelsohn (2006): http://www.cimec.org.ar/ojs/index.php/mc/article/view/486/464
582 .. _Cauchy stress tensor: https://en.wikipedia.org/wiki/Cauchy_stress_tensor
583 .. _`Chorin's projection method`: https://en.wikipedia.org/wiki/Projection_method_(fluid_dynamics)#Chorin.27s_projection_method
584 .. _`Chorin (1968)`: http://www.ams.org/journals/mcom/1968-22-104/S0025-5718-1968-0242392-2/S0025-5718-1968-0242392-2.pdf
585 .. _split: http://www.wolframalpha.com/input/?i=div(p+v)
586 .. _Poisson equation: https://en.wikipedia.org/wiki/Poisson's_equation
587 .. _`Jacobi iterative solution procedure`: http://www.rsmas.miami.edu/personal/miskandarani/Courses/MSC321/Projects/prjpoisson.pdf
588 .. _iteratively: https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
589