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