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