You do not have permission to edit this page, for the following reason:

You are not allowed to execute the action you have requested.


You can view and copy the source of this page.

x
 
1
==Lagrangian analysis of multiscale particulate flows with the particle finite element method==
2
3
'''E. Oñate'''<sup>1,2</sup>, '''M.A. Celigueta<sup>1,2</sup>''', '''S. Latorre'''<sup>1</sup>, '''G. Casas'''<sup>1,2</sup>, '''R. Rossi'''<sup>1,2</sup>, '''J. Rojek'''<sup>3</sup> 
4
5
<sup>1</sup>  Centre Internacional de Metodes Numerics a l'Enginyeria - CIMNE, Barcelona, Spain<br />
6
<sup>2</sup> Universitat Politècnica de Catalunya (UPC), Barcelona, Spain<br />
7
<sup>3</sup> IPPPT, Polish Academy of Sciences, Warsaw, Poland
8
9
==Abstract==
10
11
We present a Lagrangian numerical technique for the analysis of flows incorporating physical particles of different sizes. The numerical approach is based on the Particle Finite Element Method (PFEM) which blends concepts from particle-based techniques and the FEM. The basis of the Lagrangian formulation for particulate flows and the procedure for modelling the motion of small and large particles that are submerged in the fluid are described in detail. The numerical technique for analysis of this type of multiscale particulate flows using a stabilized mixed velocity-pressure formulation and the PFEM is also presented. Examples of application of the PFEM to several particulate flows problems are given.
12
13
'''keywords''': Lagrangian analysis, Multiscale particulate flows, Particle finite element method
14
15
==1 Introduction==
16
17
The study of fluid flows containing particles of different sizes (hereafter called ''particulate flows'') is relevant to many areas of engineering and applied sciences. In this work we are concerned with particulate flows containing small to large particles. This type of flows is typical in slurry flows originated by natural hazards such as floods, tsunamis and landslides, as well as in many processes of the bio-medical and pharmaceutical industries, in the manufacturing industry and in the oil and gas industry (i.e. cuttings transport in boreholes), among other applications <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-6'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-16'></span><span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-26'></span><span id='citeF-47'></span><span id='citeF-50'></span><span id='citeF-51'></span><span id='citeF-55'></span><span id='citeF-61'></span><span id='citeF-62'></span>[[#cite-1|[1,2,6,13,14,16,21,22,23,26,47,50,51,55,61,62]]].
18
19
Our interest in this work is the modelling and simulation of free surface particulate quasi-incompressible flows containing particles of different sizes using a particular class of Lagrangian FEM termed the Particle Finite Element Method (PFEM, www.cimne.com/pfem) [CaOnSu10,CaOnSu13,CrFrPe11,FrOnCa13,IdOnPi04,IdMaLiOn08,IdMiOn09a,IdOn10,
20
Laetal08,Oliver07,Oliver09,
21
Onetal04b,OnCeId06a,Onetal06c,Onetal08,Onetal10,Onetal11,OnCa13,OnFrCa14a,OnFrCa14b]. The PFEM treats the mesh nodes in the fluid and solid domains as particles which can freely move and even separate from the main fluid domain. A mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM.
22
23
In Lagrangian analysis procedures (such as the PFEM) the motion of fluid particles is tracked during the transient solution. Hence, the convective terms vanish in the momentum equations and no numerical stabilization is needed. Another source of instability, however, remains in the numerical solution of Lagrangian flows, that due to the treatment of the incompressibility constraint which requires using a stabilized numerical method.
24
25
In this work we use a stabilized Lagrangian formulation that has excellent mass preservation features. The success of the formulation relies on the consistent derivation of a residual-based stabilized expression of the mass balance equation using the Finite Calculus (FIC) method <span id='citeF-29'></span><span id='citeF-30'></span><span id='citeF-31'></span><span id='citeF-32'></span><span id='citeF-33'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-29|[29,30,31,32,33,37,38,39]]].
26
27
The lay-out of the paper is the following. In the next section we present the basic equations for conservation of linear momentum and mass for a quasi-incompressible particulate fluid in a Lagrangian framework. The treatment of the different force terms for micro, macro and large particles are explained. Next we derive the stabilized FIC form of the mass balance equation. Then, the finite element discretization using simplicial element with equal order approximation for the velocity and the pressure is presented and the relevant matrices and vectors of the discretized problem are given.  Details of the implicit transient solution of the Lagrangian FEM equations for a particulate flow using a Newton-Raphson type iterative scheme are presented. The basic steps of the PFEM for solving free-surface particulate flow problems are described.
28
29
The efficiency and accuracy of the PFEM  for analysis of particulate flows are verified by solving a set of  free surface and confined fluid flow problems incorporating particles of small and large sizes in two (2D) and three (3D) dimensions. The problems include the study of soil erosion, landslide situations, tsunami and flood flows, soil dredging problems and particle filling of fluid containers, among others. The excellent performance of the numerical method proposed for analysis of particulate flows is highlighted.
30
31
==2 Modelling of micro, macro and large particles==
32
33
Figure [[#img-1|1]] shows a domain containing a fluid and particles of different sizes. Particles will be termed ''microscopic'', ''macroscopic'' and ''large'' in terms of their dimensions. Microscopic and macroscopic particles will be assumed to have a cylindrical (in 2D) or spherical (in 3D) shape.  These particles are modelled as rigid objects that undergo interaction forces computed in terms of the relative distance between particles (for microscropic particles) or via the physical contact   between a particle and its neighbors (for macroscopic particles), as in the standard discrete element method (DEM) <span id='citeF-2'></span><span id='citeF-16'></span><span id='citeF-34'></span>[[#cite-2|[2,16,34]]].
34
35
<div id='img-1'></div>
36
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 70%;max-width: 100%;"
37
|-
38
|[[Image:draft_Samper_928603926-Imagen1.png|420px|Microscopic, macroscopic and large particles within a fluid domain]]
39
|- style="text-align: center; font-size: 75%;"
40
| colspan="1" | '''Figure 1:''' Microscopic, macroscopic and large particles within a fluid domain
41
|}
42
43
<div id='img-2'></div>
44
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
45
|-
46
|[[Image:draft_Samper_928603926-Imagen3.png|400px|]]
47
|[[Image:draft_Samper_928603926-Imagen4.png|400px|(a) Particles of different sizes surrounding the nodes in a FEM mesh. (b) Representative volume for a node (in shadowed  darker colour)]]
48
|- style="text-align: center; font-size: 75%;"
49
| colspan="2" | '''Figure 2:''' (a) Particles of different sizes surrounding the nodes in a FEM mesh. (b) Representative volume for a node (in shadowed  darker colour)
50
|}
51
52
In this work microscopic and macroscopic particles are assumed to be spherical and submerged in the fluid (Figure [[#img-2|2]]). Fluid-to-particle forces are transferred to the particles via appropriate drag and buoyancy functions. Particle-to-fluid forces have equal magnitude and opposite direction as the fluid-to-particle ones and are transferred to the fluid points as an additional body force vector in the momentum equation (Figure [[#img-3|3]]). These equations, as well as the mass balance equations account for the percentage of particles in the fluid, similarly at it is done in standard immersed approaches for particulate flows <span id='citeF-53'></span><span id='citeF-54'></span><span id='citeF-56'></span>[[#cite-53|[53,54,56]]].
53
54
<div id='img-3'></div>
55
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;max-width: 100%;"
56
|-
57
|[[Image:draft_Samper_928603926-Imagen5.png|540px|Immersed approach for treating the motion of physical particles in a fluid <span id='citeF-61'></span>[[#cite-61|[61]]]]]
58
|- style="text-align: center; font-size: 75%;"
59
| colspan="1" | '''Figure 3:''' Immersed approach for treating the motion of physical particles in a fluid <span id='citeF-61'></span>[[#cite-61|[61]]]
60
|}
61
62
Large particles, on the other hand, can have any arbitrary shape and they can be treated as rigid or deformable bodies. In the later case, they are discretized with the standard FEM. The forces between the fluid and the particles and viceversa are computed via fluid-structure interaction (FSI) procedures <span id='citeF-31'></span><span id='citeF-60'></span>[[#cite-31|[31,60]]].
63
64
The following sections describe the basic governing equations for a Lagrangian particulate fluid and the computation of the forces for microscopic, macroscopic and large particles.
65
66
==3 Basic governing equations for a Lagrangian particulate fluid <span id='citeF-1'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-61'></span>[[#cite-1|[1,22,23,61]]]==
67
68
==='''Conservation of linear momentum'''===
69
70
<span id="eq-1"></span>
71
{| class="formulaSCP" style="width: 100%; text-align: left;" 
72
|-
73
| 
74
{| style="text-align: left; margin:auto;" 
75
|-
76
| style="text-align: center;" | <math>\rho _f \frac{Dv_i}{Dt}=  {\partial \sigma _{ij} \over \partial x_j} +b_i + \frac{1}{n_f}  f_i^{pf}\quad , \quad i,j=1,\cdots ,n_s \quad  \hbox{in } V  </math>
77
|}
78
| style="width: 5px;text-align: right;" | (1)
79
|}
80
81
In <math display="inline">V</math> is the analysis domain, <math display="inline">n_s</math> is the number of space dimensions (<math display="inline">n_s=3</math> for 3D problems), <math display="inline">\rho _f </math> is the density of the fluid, <math display="inline">v_i</math> and <math display="inline">b_i</math> are the velocity and body force components along the <math display="inline">i</math>th cartesian axis, respectively, <math display="inline">\sigma _{ij}</math> are the fluid Cauchy stresses, <math display="inline">f_i^{pf}</math> are averaged particle-to-fluid interaction forces for which closure relations must be provided and <math display="inline">n_f</math> is the fluid volume fraction defined for each node <math display="inline">j</math> as
82
83
<span id="eq-2"></span>
84
{| class="formulaSCP" style="width: 100%; text-align: left;" 
85
|-
86
| 
87
{| style="text-align: left; margin:auto;" 
88
|-
89
| style="text-align: center;" | <math>n_{f_j} = 1- \frac{1}{V_j}\sum \limits _{i=1}^{n_j}V_j^i  </math>
90
|}
91
| style="width: 5px;text-align: right;" | (2)
92
|}
93
94
where <math display="inline">V_j</math> is the volume of the representative domain associated to a fluid node <math display="inline">j</math>, <math display="inline">V_j^i</math> is the volume of the <math display="inline">i</math>th particle belonging to <math display="inline">V_j</math> and <math display="inline">n_j</math> is the number of particles contained in <math display="inline">V_j</math>. Note that <math display="inline">n_{f_j}=1 </math> for a representative fluid domain containing no particles (Figure [[#img-2|2]]).
95
96
The fluid volume fraction <math display="inline">n_f</math> in Eq.([[#eq-1|1]]) is a continuous function that is interpolated from the nodal values in the finite element fashion <span id='citeF-22'></span>[[#cite-22|[22]]].
97
98
Summation of terms with repeated indices  is assumed in Eq.([[#eq-1|1]]) and in  the following, unless otherwise specified.
99
100
'''Remark 1.''' The term <math display="inline">\frac{Dv_i}{Dt}</math> in Eq.([[#eq-1|1]]) is the ''material derivative'' of the velocity <math display="inline">v_i</math>. This term is typically computed in a Lagrangian framework as
101
102
<span id="eq-3"></span>
103
{| class="formulaSCP" style="width: 100%; text-align: left;" 
104
|-
105
| 
106
{| style="text-align: left; margin:auto;" 
107
|-
108
| style="text-align: center;" | <math>\frac{Dv_i}{Dt} = \frac{{}^{n+1} v_i - {}^n v_i }{\Delta t}     </math>
109
|}
110
| style="width: 5px;text-align: right;" | (3)
111
|}
112
113
with
114
115
<span id="eq-4"></span>
116
{| class="formulaSCP" style="width: 100%; text-align: left;" 
117
|-
118
| 
119
{| style="text-align: left; margin:auto;" 
120
|-
121
| style="text-align: center;" | <math>{}^{n+1} v_i:=v_i ({}^{n+1} {x},{}^{n+1} t)\quad ,\quad {}^n v_i:=v_i ({}^n{x},{}^n t)       </math>
122
|}
123
| style="width: 5px;text-align: right;" | (4)
124
|}
125
126
where <math display="inline">{}^n v_i({}^n{x},{}^nt)</math> is the velocity of the material point that has the position <math display="inline">{}^n{x}</math> at time <math display="inline">t={}^nt</math>, where <math display="inline">{x}=[x_1,x_2,x_3]^T</math> is the coordinates vector of a point  in a fixed Cartesian  system. Note that the convective term, typical of Eulerian formulations, does not appear in the definition of the material derivative <span id='citeF-3'></span><span id='citeF-9'></span><span id='citeF-60'></span>[[#cite-3|[3,9,60]]].
127
128
==='''Constitutive equations'''===
129
130
The Cauchy stresses <math display="inline">\sigma _{ij}</math>  in the fluid are split in the deviatoric (<math display="inline">s_{ij}</math>) and pressure (<math display="inline">p</math>) components as
131
132
<span id="eq-5"></span>
133
{| class="formulaSCP" style="width: 100%; text-align: left;" 
134
|-
135
| 
136
{| style="text-align: left; margin:auto;" 
137
|-
138
| style="text-align: center;" | <math>\sigma _{ij} = s_{ij}+ p \delta _{ij} </math>
139
|}
140
| style="width: 5px;text-align: right;" | (5)
141
|}
142
143
where <math display="inline">\delta _{ij}</math> is the Kronecker delta. In this work the pressure is assumed to be positive for a tension state.
144
145
The relationship between the deviatoric stresses and the strain rates has the standard form for a Newtonian fluid <span id='citeF-9'></span>[[#cite-9|[9]]],
146
147
<span id="eq-6"></span>
148
{| class="formulaSCP" style="width: 100%; text-align: left;" 
149
|-
150
| 
151
{| style="text-align: left; margin:auto;" 
152
|-
153
| style="text-align: center;" | <math>s_{ij} = 2\mu \varepsilon '_{ij} \quad \hbox{with } \varepsilon '_{ij}= \varepsilon _{ij} - {1 \over 3}\varepsilon _{v} \delta _{ij} \quad \hbox{and}\quad \varepsilon _{v} = \varepsilon _{ii} </math>
154
|}
155
| style="width: 5px;text-align: right;" | (6)
156
|}
157
158
In Eq.([[#eq-6|6]]) <math display="inline">\mu </math> is the viscosity, <math display="inline">\varepsilon '_{ij}</math> and <math display="inline">\varepsilon _{v}</math> are the deviatoric and volumetric strain rates, respectively. The strain rates <math display="inline">\varepsilon _{ij}</math> are related to the velocities by
159
160
<span id="eq-7"></span>
161
{| class="formulaSCP" style="width: 100%; text-align: left;" 
162
|-
163
| 
164
{| style="text-align: left; margin:auto;" 
165
|-
166
| style="text-align: center;" | <math>\varepsilon _{ij} = {1 \over 2}\left({\partial v_i \over \partial x_j} +  {\partial v_j \over \partial x_i}  \right) </math>
167
|}
168
| style="width: 5px;text-align: right;" | (7)
169
|}
170
171
==='''Mass conservation equation'''===
172
173
The  mass conservation equation for a particulate flow is written as
174
175
<span id="eq-8"></span>
176
{| class="formulaSCP" style="width: 100%; text-align: left;" 
177
|-
178
| 
179
{| style="text-align: left; margin:auto;" 
180
|-
181
| style="text-align: center;" | <math>r_v =0  </math>
182
|}
183
| style="width: 5px;text-align: right;" | (8)
184
|}
185
186
with
187
188
<span id="eq-9"></span>
189
{| class="formulaSCP" style="width: 100%; text-align: left;" 
190
|-
191
| 
192
{| style="text-align: left; margin:auto;" 
193
|-
194
| style="text-align: center;" | <math>r_v:=\frac{D (n_f \rho _f)}{Dt}+n_f \rho _f \varepsilon _{v} </math>
195
|}
196
| style="width: 5px;text-align: right;" | (9)
197
|}
198
199
Expanding the material derivative and dividing Eq.([[#eq-8|8]]) by <math display="inline">n_f</math> the expression of <math display="inline">r_v</math> can be rewritten as
200
201
<span id="eq-10"></span>
202
{| class="formulaSCP" style="width: 100%; text-align: left;" 
203
|-
204
| 
205
{| style="text-align: left; margin:auto;" 
206
|-
207
| style="text-align: center;" | <math>r_v:= \frac{1}{\kappa } \frac{Dp}{Dt}+ \frac{1}{n_f} \frac{Dn_f}{Dt} +\varepsilon _{v}  </math>
208
|}
209
| style="width: 5px;text-align: right;" | (10)
210
|}
211
212
where <math display="inline">\kappa = \rho _f c^2</math> and <math display="inline">c= - {\partial p \over \partial \rho }</math> is the speed of sound in the fluid.
213
214
'''Remark 2.''' For <math display="inline">n_f =1</math> (i.e. no particles are contained in the fluid) the standard momentum and mass conservation equations for the fluid are recovered.
215
216
==='''Boundary conditions'''===
217
218
The boundary conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann (<math display="inline">\Gamma _t</math>) boundaries with <math display="inline">\Gamma  = \Gamma _v \cup \Gamma _t</math> are
219
220
<span id="eq-11"></span>
221
{| class="formulaSCP" style="width: 100%; text-align: left;" 
222
|-
223
| 
224
{| style="text-align: left; margin:auto;" 
225
|-
226
| style="text-align: center;" | <math>v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v </math>
227
|}
228
| style="width: 5px;text-align: right;" | (11)
229
|}
230
231
<span id="eq-12"></span>
232
{| class="formulaSCP" style="width: 100%; text-align: left;" 
233
|-
234
| 
235
{| style="text-align: left; margin:auto;" 
236
|-
237
| style="text-align: center;" | <math>\sigma _{ij}n_j -t_i^p =0 \qquad \hbox{on }\Gamma _t \quad i,j=1,\cdots ,n_s  </math>
238
|}
239
| style="width: 5px;text-align: right;" | (12)
240
|}
241
242
where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math> are the prescribed velocities and prescribed tractions on the <math display="inline">\Gamma _v</math> and <math display="inline">\Gamma _t</math> boundaries, respectively and <math display="inline">n_j</math> are the components of the unit normal vector to the boundary <span id='citeF-3'></span><span id='citeF-9'></span><span id='citeF-60'></span>[[#cite-3|[3,9,60]]].
243
244
At a free surface the Neumann boundary conditions typically apply.
245
246
==4 Motion of microscopic and macroscopic particles==
247
248
As mentioned early, microscopic and macroscopic particles are assumed to be rigid spheres (in 3D). Their motion follows the standard law for Lagrangian particles. For the <math display="inline">i</math>th spherical particle
249
250
<span id="eq-13"></span>
251
{| class="formulaSCP" style="width: 100%; text-align: left;" 
252
|-
253
| 
254
{| style="text-align: left; margin:auto;" 
255
|-
256
| style="text-align: center;" | <math>m_i \dot {u}_i = {F}_i\quad ,\quad J_i \dot {w}_i = {T}_i  </math>
257
|}
258
| style="width: 5px;text-align: right;" | (13)
259
|}
260
261
where <math display="inline">{u}_i</math> and <math display="inline">{w}_i</math> are the velocity and rotation vector of the center of gravity of the particle, <math display="inline">m_i</math> and <math display="inline">J_i </math> are the mass and rotational inertia of the particle, respectively and <math display="inline">{F}_i</math> and <math display="inline">{T}_i</math> are the vectors containing the forces and torques acting at the gravity center of the particle <span id='citeF-34'></span>[[#cite-34|[34]]].
262
263
The forces <math display="inline">{F}_i</math> acting on the <math display="inline">i</math>th particle are computed as
264
265
<span id="eq-14"></span>
266
{| class="formulaSCP" style="width: 100%; text-align: left;" 
267
|-
268
| 
269
{| style="text-align: left; margin:auto;" 
270
|-
271
| style="text-align: center;" | <math>{F}_i = {F}_i^w +  {F}_{i}^c + {F}_i^{fp}  </math>
272
|}
273
| style="width: 5px;text-align: right;" | (14)
274
|}
275
276
<math display="inline">{F}_i^w</math>, <math display="inline">{F}_{i}^c</math> and <math display="inline">{F}_i^{fp}</math> are the forces on the particle due to self-weight, contact interactions and fluid effects. These forces are computed as follows.
277
278
==='''Self-weight forces'''===
279
280
<span id="eq-15"></span>
281
{| class="formulaSCP" style="width: 100%; text-align: left;" 
282
|-
283
| 
284
{| style="text-align: left; margin:auto;" 
285
|-
286
| style="text-align: center;" | <math>{F}_i^w =-\rho _i \Omega _i {g}  </math>
287
|}
288
| style="width: 5px;text-align: right;" | (15)
289
|}
290
291
where <math display="inline">\rho _i</math> and <math display="inline">\Omega _i</math> are the density and the volume of the <math display="inline">i</math>th particle, respectively and <math display="inline">{g}</math> is the gravity acceleration vector.
292
293
==='''Contact forces'''===
294
295
<span id="eq-16"></span>
296
{| class="formulaSCP" style="width: 100%; text-align: left;" 
297
|-
298
| 
299
{| style="text-align: left; margin:auto;" 
300
|-
301
| style="text-align: center;" | <math>{F}_{i}^c = \sum \limits _{j=1}^{n_i} {F}_{ij}^c   </math>
302
|}
303
| style="width: 5px;text-align: right;" | (16)
304
|}
305
306
where <math display="inline">n_i</math> is the number of contact interfaces for the <math display="inline">i</math>th particle.
307
308
<span id="eq-17"></span>
309
{| class="formulaSCP" style="width: 100%; text-align: left;" 
310
|-
311
| 
312
{| style="text-align: left; margin:auto;" 
313
|-
314
| style="text-align: center;" | <math>{F}_{ij}^c = {F}_n^{ij} + {F}_s^{ij} = { F}_n^{ij}n_i + {F}_s^{ij}  </math>
315
|}
316
| style="width: 5px;text-align: right;" | (17)
317
|}
318
319
where <math display="inline">{F}_n^{ij}</math> and <math display="inline">{F}_s^{ij}</math> are the normal and tangential forces acting at the <math display="inline">i</math>th interface connecting particles <math display="inline">i</math> and <math display="inline">j</math> (Figure [[#img-4|4]]). These forces are computed in terms of the relative motion of the interacting particles as in the standard DEM <span id='citeF-2'></span><span id='citeF-16'></span><span id='citeF-34'></span>[[#cite-2|[2,16,34]]].
320
321
<div id='img-4'></div>
322
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
323
|-
324
|[[Image:draft_Samper_928603926-contact-particles.png|400px|]]
325
|[[Image:draft_Samper_928603926-contact.png|400px|]]
326
|-
327
|[[Image:draft_Samper_928603926-Fig3.png|400px|]]
328
|[[Image:draft_Samper_928603926-Fig5.png|400px|Interacting forces between microscopic (a) and macroscopic (b) particles]]
329
|- style="text-align: center; font-size: 75%;"
330
| colspan="2" | '''Figure 4:''' Interacting forces between microscopic (a) and macroscopic (b) particles
331
|}
332
333
For microscopic particles the tangential forces <math display="inline">{F}_s^{ij}</math> are neglected in Eq.([[#eq-17|17]]).
334
335
<br/>
336
337
'''Fluid-to-particle forces:''' <math display="inline">{F}_i^{fp} = {F}_i^d + {F}_i^b</math>, where <math display="inline">{F}_i^b</math> and <math display="inline">{F}_i^d</math> are, respectively, the buoyancy and drag forces on the <math display="inline">i</math>th particle. These forces are computed as:
338
339
==='''Buoyancy forces'''===
340
341
<span id="eq-18"></span>
342
{| class="formulaSCP" style="width: 100%; text-align: left;" 
343
|-
344
| 
345
{| style="text-align: left; margin:auto;" 
346
|-
347
| style="text-align: center;" | <math>{F}_i^b = {\Omega }_i {\boldsymbol \nabla } p  </math>
348
|}
349
| style="width: 5px;text-align: right;" | (18)
350
|}
351
352
==='''Drag forces'''===
353
354
{| class="formulaSCP" style="width: 100%; text-align: left;" 
355
|-
356
| 
357
{| style="text-align: left; margin:auto;" 
358
|-
359
| style="text-align: center;" | <math>{F}_i^d = {f}_i^d n_f^{-(\chi +1)}</math>
360
|}
361
|}
362
363
where <math display="inline">\chi = \chi (Re)</math> is a coefficient that depends on the local Reynolds number for the particle <math display="inline">Re</math> <span id='citeF-1'></span><span id='citeF-6'></span><span id='citeF-15'></span><span id='citeF-22'></span><span id='citeF-23'></span>[[#cite-1|[1,6,15,22,23]]] and
364
365
<span id="eq-19"></span>
366
{| class="formulaSCP" style="width: 100%; text-align: left;" 
367
|-
368
| 
369
{| style="text-align: left; margin:auto;" 
370
|-
371
| style="text-align: center;" | <math>{f}_i^d = \frac{1}{2} \rho _f A_i C_d \Vert {v}_{f_i} - {v}_i\Vert ({v}_{f_i} - {v}_i)  </math>
372
|}
373
| style="width: 5px;text-align: right;" | (19)
374
|}
375
376
In Eq.([[#eq-19|19]]) <math display="inline">{v}_{fi}</math> and <math display="inline">{v}_i</math> are respectively the velocity of the fluid and of the particle center, <math display="inline">A_i</math> is the area of the particle surface with radius <math display="inline">r_i</math> (<math display="inline">2\pi r_i</math> or <math display="inline">4\pi r_i^2</math> for a circle or a sphere, respectively) and <math display="inline">C_d</math> is a drag coefficient that depends on the particle geometry and the rugosity of its surface <span id='citeF-22'></span><span id='citeF-23'></span>[[#cite-22|[22,23]]].
377
378
The force term <math display="inline"> {f}^{pf}_i</math> in the r.h.s. of the momentum equations (Eq.([[#eq-1|1]])) is computed for each particle (in vector form) as <math display="inline"> {f}^{pf}=-{f}^{fp}</math> with vector <math display="inline"> {f}^{fp}</math> computed ''at each node'' in the fluid mesh         from the drag forces <math display="inline">{F}_i^{d}</math> as
379
380
<span id="eq-20"></span>
381
{| class="formulaSCP" style="width: 100%; text-align: left;" 
382
|-
383
| 
384
{| style="text-align: left; margin:auto;" 
385
|-
386
| style="text-align: center;" | <math>{f}^{fp}_j = \frac{1}{V_j} \sum \limits _{i=1}^{n_j} {F}_i^d~~~,~~ j=1,N  </math>
387
|}
388
| style="width: 5px;text-align: right;" | (20)
389
|}
390
391
A continuum distribution of <math display="inline">{f}^{fp}</math> is obtained by interpolating its nodal values over each element in the FEM fashion.
392
393
We note that the forces on the particles due to lift effects have been neglected in the present analysis. These forces can be accounted for as explained in <span id='citeF-22'></span>[[#cite-22|[22]]].
394
395
==5 Motion of large particles==
396
397
As mentioned earlier, large particles may be considered as rigid or deformable bodies. In the first case the motion follows the rules of Eq.([[#eq-13|13]]) for rigid Lagrangian particles. The contact forces at the particle surface due to the adjacent interacting particles are computed using a frictional contact interface layer between particles as in the standard PFEM (Section 10.2).
398
399
The fluid forces on the particles are computing by integrating the tangential (viscous) and normal (pressure) forces at the edges of the fluid elements surrounding the particles.
400
401
Large deformable particles, on the other hand, behave as deformable bodies immersed in the fluid which are discretized via the standard FEM. Their motion can be followed using a staggered FSI scheme, or else by treating the deformable bodies and the fluid as a single continuum with different material properties. Details of this unified treatment of the interaction between fluids and deformable solids can be found in <span id='citeF-12'></span><span id='citeF-18'></span><span id='citeF-46'></span>[[#cite-12|[12,18,46]]].
402
403
==6 Stabilized FIC form of the mass balance equation==
404
405
The modelling of incompressible fluids with a mixed finite element method using an equal order interpolation for the velocities and the pressure requires introducing a stabilized formulation for the mass balance equation.
406
407
In our work we use a stabilized form of the mass balance equation obtained via the Finite Calculus (FIC) approach <span id='citeF-29'></span><span id='citeF-30'></span><span id='citeF-31'></span><span id='citeF-32'></span><span id='citeF-33'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-29|[29,30,31,32,33,37,38,39]]]. The FIC stabilized mass balance equation is written as
408
409
<span id="eq-21"></span>
410
{| class="formulaSCP" style="width: 100%; text-align: left;" 
411
|-
412
| 
413
{| style="text-align: left; margin:auto;" 
414
|-
415
| style="text-align: center;" | <math>r_v - \tau {\partial \bar r_{m_i} \over \partial x_i}=0 \quad \hbox{in } V  </math>
416
|}
417
| style="width: 5px;text-align: right;" | (21)
418
|}
419
420
where
421
422
<span id="eq-22"></span>
423
{| class="formulaSCP" style="width: 100%; text-align: left;" 
424
|-
425
| 
426
{| style="text-align: left; margin:auto;" 
427
|-
428
| style="text-align: center;" | <math>\bar r_{m_i} = {\partial \sigma _{ij} \over \partial x_j}+b_i   + \frac{1}{n_f} f_i^{pf}  </math>
429
|}
430
| style="width: 5px;text-align: right;" | (22)
431
|}
432
433
is a static momentum term and <math display="inline">\tau </math> is a stabilization parameter computed as
434
435
<span id="eq-23"></span>
436
{| class="formulaSCP" style="width: 100%; text-align: left;" 
437
|-
438
| 
439
{| style="text-align: left; margin:auto;" 
440
|-
441
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{h^2}+\frac{2\rho _f}{\delta }   \right)^{-1}  </math>
442
|}
443
| style="width: 5px;text-align: right;" | (23)
444
|}
445
446
where <math display="inline">h</math> is a characteristic distance of each finite element and <math display="inline">\delta </math> is a time parameter.
447
448
The derivation of Eq.([[#eq-21|21]]) for an homogeneous quasi-incompressible fluid is presented in <span id='citeF-45'></span>[[#cite-45|[45]]].
449
450
The stabilization parameter <math display="inline">\tau </math>  is computed in practice for each element <math display="inline">e</math> using <math display="inline">h=l^e</math> and <math display="inline">\delta = \Delta t</math> as
451
452
<span id="eq-24"></span>
453
{| class="formulaSCP" style="width: 100%; text-align: left;" 
454
|-
455
| 
456
{| style="text-align: left; margin:auto;" 
457
|-
458
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{(l^e)^2} + \frac{2\rho }{\Delta t}  \right)^{-1} </math>
459
|}
460
| style="width: 5px;text-align: right;" | (24)
461
|}
462
463
where <math display="inline">\Delta t</math> is the time step used for the transient solution and <math display="inline">l^e</math> is a characteristic element length computed as <math display="inline">l^e = 2(V^e)^{1/n_s}</math> where <math display="inline">V^e</math> is the element area (for 3-noded triangles) or volume (for 4-noded tetrahedra). For fluids with heterogeneous material properties the values of <math display="inline">\mu </math> and <math display="inline">\rho </math> in Eq.([[#eq-24|24]]) are computed at the element center.
464
465
==7 Variational equations for the fluid==
466
467
The variational form of the momentum and mass balance equations is obtained via the standard weighted residual approach <span id='citeF-9'></span><span id='citeF-60'></span>[[#cite-9|[9,60]]]. The resulting integral expressions after integration by parts and some algebra are:
468
469
==='''Momentum equations'''===
470
471
<span id="eq-25"></span>
472
{| class="formulaSCP" style="width: 100%; text-align: left;" 
473
|-
474
| 
475
{| style="text-align: left; margin:auto;" 
476
|-
477
| style="text-align: center;" | <math>\int _V \! w_i \rho \frac{Dv_i}{Dt} dV + \!\int _V \left[\delta \varepsilon _{ij} 2\mu \varepsilon '_{ij} + \delta \varepsilon _{v}p\right]dV - \!\int _V w_i \left(b_i   + \frac{1}{n_f} f_i^{pf}\right)dV - \!\int _{\Gamma _t} w_i t_i^p d\Gamma =0  </math>
478
|}
479
| style="width: 5px;text-align: right;" | (25)
480
|}
481
482
==='''Mass balance equation'''===
483
484
<span id="eq-26"></span>
485
{| class="formulaSCP" style="width: 100%; text-align: left;" 
486
|-
487
| 
488
{| style="text-align: left; margin:auto;" 
489
|-
490
| style="text-align: center;" | <math>\begin{array}{r}\displaystyle \int _V \frac{q}{\kappa } \frac{Dp}{Dt}dV - \int _V q \left(\frac{1}{n_f} \frac{Dn_f}{Dt} +\varepsilon _{v}\right)dV + \int _V \tau {\partial q \over \partial x_i} \left({\partial  \over \partial x_i} (2\mu \varepsilon _{ij})+ {\partial p \over \partial x_i}+b_i\right)dV\\ \displaystyle - \int _{\Gamma _t} q \tau \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} (2\mu {\partial v_n \over \partial n} + p -t_n)\right]d\Gamma =0 \end{array}   </math>
491
|}
492
| style="width: 5px;text-align: right;" | (26)
493
|}
494
495
The derivation of Eqs.([[#eq-25|25]]) and ([[#eq-26|26]]) for homogeneous fluids can be found in <span id='citeF-45'></span>[[#cite-45|[45]]].
496
497
==8 FEM discretization==
498
499
We discretize  the analysis domain containing <math display="inline">N_p</math> microscopic and macroscopic particles and a number of large particles into finite elements with <math display="inline">n</math> nodes in the standard manner leading to a mesh with a total number of <math display="inline">N_e</math> elements and <math display="inline">N</math> nodes. In our work we will choose simple 3-noded linear triangles (<math display="inline">n=3</math>) for 2D problems and 4-noded tetrahedra (<math display="inline">n=4</math>) for 3D problems with local linear shape functions <math display="inline">N_i^e</math> defined for each node <math display="inline">i</math> of element <math display="inline">e</math> <span id='citeF-41'></span><span id='citeF-58'></span>[[#cite-41|[41,58]]]. The velocity components, the weighting functions and the pressure are interpolated over the mesh in terms of their nodal values in the same manner using the  global linear shape functions <math display="inline">N_j</math> spanning over the elements sharing node <math display="inline">j</math> (<math display="inline">j=1,N</math>)  <span id='citeF-41'></span><span id='citeF-58'></span>[[#cite-41|[41,58]]].
500
501
The finite element interpolation over the fluid domain can be written in matrix form as
502
503
<span id="eq-27"></span>
504
{| class="formulaSCP" style="width: 100%; text-align: left;" 
505
|-
506
| 
507
{| style="text-align: left; margin:auto;" 
508
|-
509
| style="text-align: center;" | <math>{v} = {N}_v\bar {v} ~,~{w} = {N}_v\bar {w} ~,~{p} = {N}_p\bar {p}~,~ q= {N}_p \bar {q}  </math>
510
|}
511
| style="width: 5px;text-align: right;" | (27)
512
|}
513
514
where
515
516
<span id="eq-28"></span>
517
{| class="formulaSCP" style="width: 100%; text-align: left;" 
518
|-
519
| 
520
{| style="text-align: left; margin:auto;" 
521
|-
522
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \bar {v} = \left\{\begin{matrix}\bar{v}^1\\\bar{v}^2\\\vdots \\ \bar{v}^N\end{matrix}  \right\}\quad \hbox{with}~ \bar{v}^i = \left\{\begin{matrix}\bar{v}^i_1\\\bar{v}^i_2\\\bar{v}^i_3\end{matrix}  \right\}~, ~ \bar {w}= \left\{\begin{matrix}\bar{w}^1\\\bar{w}^2\\\vdots \\ \bar{w}^N\end{matrix}  \right\}\quad \hbox{with}~ \bar{w}^i = \left\{\begin{matrix}\bar{w}^i_1\\\bar{w}^i_2\\\bar{w}^i_3\end{matrix}  \right\}~,~\bar {p} =\left\{\begin{matrix}\bar{p}^1\\\bar{p}^2\\\vdots \\ \bar{p}^{N}\end{matrix}  \right\}\quad \hbox{and}~ \bar {q} =\left\{\begin{matrix}\bar{q}^1\\\bar{q}^2\\\vdots \\ \bar{q}^{N}\end{matrix}  \right\}\\ \displaystyle {N}_v = [{N}_1, {N}_2,\cdots , {N}_N ]^T ~,~ {N}_p = [{N}_1, { N}_2,\cdots , {N}_N ] \end{array} </math>
523
|}
524
| style="width: 5px;text-align: right;" | (28)
525
|}
526
527
with <math display="inline">{N}_j = N_j {I}_{n_s}</math> where <math display="inline">{I}_{n_s}</math> is the <math display="inline">n_s\times n_s</math> unit matrix.
528
529
In Eq.([[#eq-28|28]]) vectors <math display="inline">\bar{v}</math>, <math display="inline">(\bar{w},\bar {q})</math> and <math display="inline">\bar {p}</math> contain the nodal velocities, the nodal weighting functions and the nodal pressures for the whole mesh, respectively and the upperindex denotes the nodal value for each vector or scalar magnitude.
530
531
Substituting the approximation ([[#eq-27|27]]) into the variational forms ([[#eq-25|25]]) and ([[#eq-26|26]]) gives the system of algebraic equations for the particulate fluid in matrix form as
532
533
<span id="eq-29"></span>
534
{| class="formulaSCP" style="width: 100%; text-align: left;" 
535
|-
536
| 
537
{| style="text-align: left; margin:auto;" 
538
|-
539
| style="text-align: center;" | <math>{M}_0 {\dot{\bar{v}}} + {K}\bar{v}+{Q}\bar {p}- {f}_v={0} </math>
540
|}
541
| style="width: 5px;text-align: right;" | (29)
542
|}
543
544
<span id="eq-30"></span>
545
{| class="formulaSCP" style="width: 100%; text-align: left;" 
546
|-
547
| 
548
{| style="text-align: left; margin:auto;" 
549
|-
550
| style="text-align: center;" | <math>\hbox{M}_1 {\dot{\bar{p}}}-{Q}^T \bar{v} + ({L}+{M}_b)\bar {p}- {f}_p={0} </math>
551
|}
552
| style="width: 5px;text-align: right;" | (30)
553
|}
554
555
The different matrices and vectors in Eqs.([[#8 FEM discretization|8]]) are shown in Box 1 for 2D problems.
556
557
'''Remark 3.''' The boundary terms of vector <math display="inline">{f}_p</math> can be incorporated in the matrices of Eq.([[#eq-30|30]]). This, however, leads to a non symmetrical set of equations. For this reason we have chosen to compute these boundary terms iteratively within the incremental solution scheme.
558
559
'''Remark 4.''' Matrix <math display="inline">{M}_b</math> in Eq.([[#eq-30|30]]) allows us to compute the pressure without the need of prescribing its value at the free surface. This eliminates the error introduced when the pressure is prescribed to zero in free boundaries, which may lead to considerable mass losses <span id='citeF-20'></span><span id='citeF-45'></span>[[#cite-20|[20,45]]].
560
561
==9 Incremental solution of the discretized equations==
562
563
Eqs.([[#8 FEM discretization|8]]) are solved in time with an implicit Newton-Raphson type iterative scheme <span id='citeF-3'></span><span id='citeF-9'></span><span id='citeF-58'></span><span id='citeF-60'></span>[[#cite-3|[3,9,58,60]]]. The basic steps within a time interval <math display="inline">[n,n+1]</math> are:<br/>
564
565
* Initialize variables: <math display="inline">({}^{n+1}\hbox{x}^1,{}^{n+1}\bar {v}^1,{}^{n+1}\bar{p}^1,{}^{n+1}n_f^i,   {}^{n+1}\bar{r}^1_m)\equiv  \left\{{}^{n}{x},{}^{n}\bar{v},{}^{n}\bar{p},{}^{n} n_f ,{}^{n}\bar{r}_m \right\}</math>.
566
* Iteration loop: <math display="inline">k=1,\cdots , NITER</math>.<p> For each iteration. 
567
568
</p>
569
570
<br/>
571
572
===''Step 1. Compute the nodal velocity increments <math>\Delta \bar {v}</math>''===
573
574
From Eq.([[#8 FEM discretization|8]]a), we deduce
575
576
<span id="eq-31"></span>
577
{| class="formulaSCP" style="width: 100%; text-align: left;" 
578
|-
579
| 
580
{| style="text-align: left; margin:auto;" 
581
|-
582
| style="text-align: center;" | <math>{}^{n+1}{H}_v^i \Delta \bar {v} = - {}^{n+1}\bar{\hbox{r}}_m^k \rightarrow  \Delta \bar{v} </math>
583
|}
584
| style="width: 5px;text-align: right;" | (31)
585
|}
586
587
with the momentum residual <math display="inline">\bar{r}_m</math> and the iteration matrix <math display="inline">{H}_v</math>  given by
588
589
<span id="eq-32"></span>
590
{| class="formulaSCP" style="width: 100%; text-align: left;" 
591
|-
592
| 
593
{| style="text-align: left; margin:auto;" 
594
|-
595
| style="text-align: center;" | <math>\bar{\hbox{r}}_m = \hbox{M}_0 {\dot{\bar{v}}} + \hbox{K}\bar{v}+\hbox{Q}\bar {p}-\hbox{f}_v\quad , \quad \hbox{H}_v = \frac{1}{\Delta t} \hbox{M}_0 + \hbox{K} + \hbox{K}_v  </math>
596
|}
597
| style="width: 5px;text-align: right;" | (32)
598
|}
599
600
with
601
602
<span id="eq-33"></span>
603
{| class="formulaSCP" style="width: 100%; text-align: left;" 
604
|-
605
| 
606
{| style="text-align: left; margin:auto;" 
607
|-
608
| style="text-align: center;" | <math>\hbox{K}_v^e = \int _{{}^nV^e} \hbox{B}^T \hbox{m} \theta \Delta t\kappa \hbox{m}^T \hbox{B}dV </math>
609
|}
610
| style="width: 5px;text-align: right;" | (33)
611
|}
612
613
===''Step 2. Update the  velocities''===
614
615
<span id="eq-34"></span>
616
{| class="formulaSCP" style="width: 100%; text-align: left;" 
617
|-
618
| 
619
{| style="text-align: left; margin:auto;" 
620
|-
621
| style="text-align: center;" | <math>\hbox{Fluid nodes:}~~{}^{n+1}\bar{v}^{k+1}= {}^{n+1}\bar{v}^k + \Delta \bar{v} </math>
622
|}
623
| style="width: 5px;text-align: right;" | (34)
624
|}
625
626
<span id="eq-35"></span>
627
{| class="formulaSCP" style="width: 100%; text-align: left;" 
628
|-
629
| 
630
{| style="text-align: left; margin:auto;" 
631
|-
632
| style="text-align: center;" | <math>\hbox{Rigid particles:}~~\left\{\begin{array}{l}{}^{n+1/2}\dot{u}_{j} = {}^{n-1/2}\dot{u}_{j} + {}^{n}\ddot{u}_{j}^{k+1} \Delta t\\ \dot{u}_{j} = \frac{1}{m_j}  {}^{n}{F}_{j}^{k+1} ~~,~~j=1,N_p \end{array}\right.  </math>
633
|}
634
| style="width: 5px;text-align: right;" | (35)
635
|}
636
637
===''Step 3. Compute the nodal pressures'' <math> {}^{n+1}\bar{p}^{k+1}</math>===
638
639
From Eq.([[#eq-30|30]]) we obtain
640
641
<span id="eq-36"></span>
642
{| class="formulaSCP" style="width: 100%; text-align: left;" 
643
|-
644
| 
645
{| style="text-align: left; margin:auto;" 
646
|-
647
| style="text-align: center;" | <math>{}^{n+1} \hbox{H}_p^i  {}^{n+1}\bar{p}^{k+1}= \frac{1}{\Delta t} \hbox{M}_1 {}^{n+1}\bar{p}^i + {Q}^T {}^{n+1}\bar{v}^{k+1}+ {}^{n+1}\bar{f}^{i}_p  \rightarrow  {}^{n+1}\bar{p}^{k+1}   </math>
648
|}
649
| style="width: 5px;text-align: right;" | (36)
650
|}
651
652
with
653
654
<span id="eq-37"></span>
655
{| class="formulaSCP" style="width: 100%; text-align: left;" 
656
|-
657
| 
658
{| style="text-align: left; margin:auto;" 
659
|-
660
| style="text-align: center;" | <math>\hbox{H}_p = \frac{1}{\Delta t} \hbox{M}_1 + \hbox{L} + \hbox{M}_b </math>
661
|}
662
| style="width: 5px;text-align: right;" | (37)
663
|}
664
665
===''Step 4. Update the coordinates of the fluid nodes and particles''===
666
667
<span id="eq-38"></span>
668
{| class="formulaSCP" style="width: 100%; text-align: left;" 
669
|-
670
| 
671
{| style="text-align: left; margin:auto;" 
672
|-
673
| style="text-align: center;" | <math>\hbox{Fluid nodes:}~~{}^{n+1}{x}^{k+1}_i= {}^{n+1}{x}^k_i + \frac{1}{2} ({}^{n+1}\bar{v}^{k+1}_i + {}^{n}\bar{v}_i)\Delta t \quad ,~~ i=1,N </math>
674
|}
675
| style="width: 5px;text-align: right;" | (38)
676
|}
677
678
<span id="eq-39"></span>
679
{| class="formulaSCP" style="width: 100%; text-align: left;" 
680
|-
681
| 
682
{| style="text-align: left; margin:auto;" 
683
|-
684
| style="text-align: center;" | <math>\hbox{Rigid particles:}~~\left\{\begin{array}{l}{}^{n+1}{u}_{i}^{k+1} = {}^{n}{u}_{i}^{k+1} + {}^{n+1/2}\dot{u}_{i}^{k+1}\Delta t\\[.4cm] {}^{n+1}{x}^{k+1}_i = {}^{n}{x}_i + {}^{n+1}{u}^{k+1}_i \quad ,~~ i=1,N_p \end{array}\right.  </math>
685
|}
686
| style="width: 5px;text-align: right;" | (39)
687
|}
688
689
===''Step 5. Compute the fluid volume fractions for each node <math>{}^{n+1}{n}^{k+1}_{f_i}</math> via Eq.(2)''===
690
691
===''Step 6. Compute forces and torques on particles:'' <math>{}^{n+1}{F}^{k+1}_i,{}^{n+1}{T}^{k+1}_i~ ,~i=1,N_p</math>===
692
693
===''Step 7. Compute particle-to-fluid nodes:'' <math>({}^{n+1}{f}^{pf}_i)^{k+1} = - ({}^{n+1}{f}^{fp}_i)^{k+1} ~ ,~ i=1,N</math> with <math>{f}^{fp}_i</math> computed by Eq.(18)===
694
695
===''Step 8. Check convergence''===
696
697
Verify the following conditions:
698
699
<span id="eq-40"></span>
700
{| class="formulaSCP" style="width: 100%; text-align: left;" 
701
|-
702
| 
703
{| style="text-align: left; margin:auto;" 
704
|-
705
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \Vert {}^{n+1}\bar {v}^{k+1}- {}^{n+1}\bar{v}^{k}\Vert \le e_v \Vert {}^{n}\bar {v}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {p}^{k+1}-{}^{n+1}\bar {p}^k\Vert \le e_p  \Vert {}^{n}\bar{p}\Vert  \end{array}   </math>
706
|}
707
| style="width: 5px;text-align: right;" | (40)
708
|}
709
710
where <math display="inline">e_v</math> and <math display="inline">e_p</math> are prescribed error norms for the nodal velocities and the nodal pressures, respectively. In the examples solved in this work we have set <math display="inline">e_v = e_p=10^{-3}</math>.
711
712
If both conditions ([[#eq-40|40]]) are satisfied then make <math display="inline">n \leftarrow n+1 </math> and proceed to the next  time step.
713
714
Otherwise, make the iteration counter <math display="inline">k \leftarrow k+1 </math> and repeat Steps 1&#8211;8.
715
716
'''Remark 5.''' In Eqs.([[#''Step 1. Compute the nodal velocity increments <math>\Delta \bar {v}</math>''|9]])&#8211;([[#eq-40|40]]) <math display="inline">{}^{n+1}(\cdot )</math> denotes the values of a matrix or a vector computed using the nodal unknowns at time <math display="inline">n+1</math>. In our work the derivatives and integrals in the iteration matrices <math display="inline">{H}_v</math> and <math display="inline">{H}_p</math> and the residual vector <math display="inline">\bar{r}_m</math>  are computed on the ''discretized geometry at time'' <math display="inline">n</math> (i.e. <math display="inline">V^e = {}^{n}V^e</math>) while the nodal force vectors <math display="inline">{f}_v</math> and <math display="inline">{f}_p</math> are computed on the current configuration at time <math display="inline">n+1</math>. This is equivalent to using an ''updated Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-12'></span><span id='citeF-44'></span><span id='citeF-59'></span>[[#cite-3|[3,12,44,59]]].
717
718
'''Remark 6.''' The tangent “bulk” stiffness matrix <math display="inline">\hbox{K}_v</math> in the iteration matrix <math display="inline">\hbox{H}_v</math> of Eq.([[#eq-32|32]]) accounts for the changes of the pressure due to the velocity. Including matrix <math display="inline">\hbox{K}_v</math> in <math display="inline">\hbox{H}_v</math> has proven to be essential for the fast convergence, mass preservation and overall accuracy of the iterative solution <span id='citeF-11'></span><span id='citeF-45'></span>[[#cite-11|[11,45]]].
719
720
'''Remark 7.'''  The parameter <math display="inline">\theta </math> in <math display="inline">{K}_v</math> (<math display="inline">0<   \theta \le 1</math>)  has the role of preventing the ill-conditioning of the iteration matrix <math display="inline">{H}_v</math> for very large values of the speed of sound in the fluid that lead to a dominant role of the terms of the tangent bulk stiffness matrix <math display="inline">{K}_v</math>. An adequate selection of <math display="inline">\theta </math> also improves the overall accuracy of the numerical solution and the preservation of mass for large time steps. Details of the derivation of Eq.(28c) can be found in <span id='citeF-45'></span>[[#cite-45|[45]]].
721
722
'''Remark 8.''' The iteration matrix <math display="inline">\hbox{H}_v</math> in Eq.([[#eq-31|31]]) is an ''approximation of the exact tangent matrix'' in the updated Lagrangian formulation for a quasi/fully incompressible fluid <span id='citeF-44'></span>[[#cite-44|[44]]]. The simplified form of <math display="inline">\hbox{H}_v</math>  used in this work has yielded very good results  with convergence achieved for the nodal velocities and pressure in 3&#8211;4 iterations in all the problems analyzed.
723
724
'''Remark 9.''' The time step within a time interval <math display="inline">[n,n+1]</math> has been chosen as <math display="inline">\Delta t =\min \left(\frac{^n l_{\min }^e}{\vert {}^n{v}\vert _{\max }},\Delta t_b\right)</math> where <math display="inline">^n l_{\min }^e</math> is the minimum characteristic distance  of all elements in the mesh, with <math display="inline">l^e</math> computed as explained in Section 6, <math display="inline">\vert ^n{v}\vert _{\max }</math> is the maximum value of the modulus of the velocity of all nodes in the mesh  and <math display="inline">\Delta t_b</math> is the critical time step of all nodes approaching a solid boundary <span id='citeF-45'></span>[[#cite-45|[45]]].
725
726
==10 About the particle finite element method (PFEM)==
727
728
===10.1 The basis of the PFEM===
729
730
Let us consider a domain <math display="inline">V</math> containing  fluid and solid subdomains. Each subdomain is characterized by a set of points, hereafter termed ''virtual particles''. The virtual particles contain all the information  for defining the geometry and the material and mechanical properties of the underlying subdomain. In the PFEM both subdomains are modelled using an ''updated'' ''Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-44'></span><span id='citeF-59'></span>[[#cite-3|[3,44,59]]].
731
732
The solution steps within a time step in the PFEM are as follows:
733
734
<ol>
735
736
<li>The starting point at each time step is the cloud of points <math display="inline">C</math> in the fluid and solid domains. For instance <math display="inline">^{n} C</math> denotes the cloud at time <math display="inline">t={}^n t </math> (Figure [[#img-5|5]]). </li>
737
738
<li>Identify the boundaries defining the analysis domain <math display="inline">^{n} V</math>, as well as the subdomains in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution, including separation and re-entering of nodes. The Alpha Shape method <span id='citeF-10'></span>[[#cite-10|[10]]] is used for the boundary definition. Clearly, the accuracy in the reconstruction of the boundaries depends on the number of points in the vicinity of each boundary and on the Alpha Shape parameter. In the problems solved in this work the Alpha Shape method has been implementation as described in <span id='citeF-17'></span><span id='citeF-35'></span>[[#cite-17|[17,35]]]. </li>
739
740
<li> Discretize the the analysis domain <math display="inline">^{n} V</math> with a finite element mesh <math display="inline">{}^{n} M.</math>We use an efficient mesh generation scheme based on an enhanced Delaunay tesselation <span id='citeF-17'></span><span id='citeF-35'></span>[[#cite-17|[17,35]]]. </li>
741
742
<li> Solve the Lagrangian equations of motion for the overall continuum using the standard FEM. Compute the state variables in at the next (updated) configuration for <math display="inline">^n t+\Delta t</math>: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. </li>
743
744
<li>  Move the mesh nodes to a new position <math display="inline">^{n+1} C</math> where ''n''+1 denotes the time <math display="inline">{}^n t +\Delta t</math>, in terms of the time increment size. </li>
745
746
<li> Go back to step 1 and repeat the solution for the next time step to obtain <math display="inline">{}^{n+2} C</math>. </li>
747
748
</ol>
749
750
<div id='img-5'></div>
751
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
752
|-
753
|[[Image:draft_Samper_928603926-Figure2a.png|600px|Sequence of steps in the PFEM to update a “cloud” of nodes representing a domain containing a fluid and a solid part (in darker shadow) from time n   (t=ⁿt)  to   time n+2 (t=ⁿt +2∆t)  ]]
754
|- style="text-align: center; font-size: 75%;"
755
| colspan="1" | '''Figure 5:''' Sequence of steps in the PFEM to update a “cloud” of nodes representing a domain containing a fluid and a solid part (in darker shadow) from time <math>n</math>   (<math>t=^n t</math>)  to   time <math>n+2</math> (<math>t=^n t +2\Delta t</math>)  
756
|}
757
758
Note that the key differences between the PFEM and the classical FEM are the remeshing technique and the identification of the domain boundary at each time step.
759
760
The CPU time required for meshing grows linearly with the number of nodes. As a general rule,  meshing consumes for 3D problems around 15% of the total CPU time per time step <span id='citeF-43'></span>[[#cite-43|[43]]].
761
762
Application of the PFEM in fluid and solid mechanics and in fluid-structure interaction problems can be found in [CaOnSu10,CaOnSu13,CrFrPe11,FrOnCa13,IdOnPi04,IdMaLiOn08,IdMiOn09a,IdOn10,
763
Laetal08,Oliver07,Oliver09,
764
Onetal04b,OnCeId06a,Onetal06c,Onetal08,Onetal10,Onetal11,OnCa13,OnFrCa14a,OnFrCa14b], as well in www.cimne.com/pfem.
765
766
===10.2 Treatment of contact conditions===
767
768
Known velocities at boundaries in the PFEM are prescribed in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. Surface tractions are applied to the Neumann part of the boundary, as usual in the FEM.
769
770
Contact between fluid particles and fixed boundaries is accounted for by adjusting the time step so that fluid nodes do not penetrate into the solid boundaries <span id='citeF-45'></span>[[#cite-45|[45]]].
771
772
The contact between two large particles (and between two bodies, in general) is treated by introducing a layer of ''contact elements'' between the two interacting particles. This contact interface layer is ''automatically created during the mesh generation step'' by prescribing a minimum distance <math display="inline">\left(h_{c} \right)</math> between two interacting particles. If the distance exceeds the minimum value <math display="inline">\left(h_{c} \right)</math> then the generated elements are treated as fluid elements. Otherwise the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacement is introduced  <span id='citeF-35'></span><span id='citeF-40'></span><span id='citeF-43'></span>[[#cite-35|[35,40,43]]] (Figure [[#img-6|6]]).
773
774
<div id='img-6'></div>
775
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
776
|-
777
|[[Image:draft_Samper_928603926-Imagen2.png|400px|]]
778
|[[Image:draft_Samper_928603926-contact-PFEM-particles.png|400px|(a) Large particles (in dark shadow) surrounded by a finite element mesh. The contact interface is shown in lighter shadow. (b) Contact interface between two objects treated as large particles and between an object and a wall]]
779
|- style="text-align: center; font-size: 75%;"
780
| colspan="2" | '''Figure 6:''' (a) Large particles (in dark shadow) surrounded by a finite element mesh. The contact interface is shown in lighter shadow. (b) Contact interface between two objects treated as large particles and between an object and a wall
781
|}
782
783
This algorithm allows us to model complex frictional contact conditions between two or more interacting bodies moving in water in an a simple manner. The algorithm has  been used to model frictional contact situations between rigid and elastic solids in structural mechanics applications, such as soil/rock excavation problems <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]]. The frictional contact algorithm described above has been extended by Oliver ''et al.'' <span id='citeF-27'></span><span id='citeF-28'></span>[[#cite-27|[27,28]]] for analysis of metal cutting and machining problems.
784
785
Figure [[#img-7|7]] shows an example of the analysis with the PFEM of a collection of large particles submerged in a tank containing water in sloshing motion.
786
787
<div id='img-7'></div>
788
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 70%;max-width: 100%;"
789
|-
790
|[[Image:draft_Samper_928603926-vlcsnap-2014-02-25-16h00m21s64.png|420px|PFEM results for the motion of large particles submerged  in a tank containing water in sloshing motion]]
791
|- style="text-align: center; font-size: 75%;"
792
| colspan="1" | '''Figure 7:''' PFEM results for the motion of large particles submerged  in a tank containing water in sloshing motion
793
|}
794
795
===10.3 Treatment of surface erosion===
796
797
Prediction of bed erosion and sediment transport in open channel flows are  important tasks in many areas of river and environmental engineering. Bed erosion can lead to instabilities of the river basin slopes. It can also undermine the foundation of bridge piles thereby favouring structural failure. Modeling of bed erosion is also relevant for predicting the evolution of surface material dragged in earth dams in overspill situations. Bed erosion is one of the main causes of environmental damage in floods.
798
799
Oñate et al. <span id='citeF-36'></span>[[#cite-36|[36]]] have proposed an extension of the PFEM to model bed erosion. The erosion model is based on the frictional work  at the bed surface originated by the shear stresses in the fluid.
800
801
The algorithm for modeling the  erosion of soil/rock particles at the fluid bed is briefly the following:
802
803
<ol>
804
805
<li>Compute at every point of the bed surface the  tangential   stress <math display="inline">\tau </math> induced by the fluid motion. </li>
806
807
<li>Compute the frictional work <math display="inline">W_f</math> originated by the tangential stress at the   bed surface. </li>
808
809
<li>The onset of erosion at a bed point occurs when <math display="inline">{}^nW_f</math> exceeds a critical   threshold value <math display="inline">W_c</math>. </li>
810
811
<li>If <math display="inline">{}^nW_f > W_c</math> at a bed node, then the node is detached from the bed   region and it is allowed to move with the fluid flow. As a consequence, the   mass of the patch of bed elements surrounding the bed node vanishes in the   bed domain and it is transferred to the adjacent   fluid node. This mass is subsequently transported with the fluid as an immersed macroscopic particle. </li>
812
813
<li>Sediment deposition can be modeled by an inverse process to that described   in the previous step. Hence, a suspended node adjacent to the bed surface with a   velocity below a threshold value is attached to the bed surface. </li>
814
815
</ol>
816
817
Figure [[#img-8|8]] shows an schematic view of the bed erosion algorithm described. Details of the algorithm can be found in <span id='citeF-36'></span>[[#cite-36|[36]]].
818
819
<div id='img-8'></div>
820
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;max-width: 100%;"
821
|-
822
|[[Image:draft_Samper_928603926-bed-erosion-2.png|540px|Modeling of bed erosion with the PFEM. The mass of the eroded domain is assigned to the fluid node k]]
823
|- style="text-align: center; font-size: 75%;"
824
| colspan="1" | '''Figure 8:''' Modeling of bed erosion with the PFEM. The mass of the eroded domain is assigned to the fluid node <math>k</math>
825
|}
826
827
==11 A nodal  algorithm for transporting  microscopic and macroscopic particles within a finite element mesh==
828
829
The fact that in the PFEM a new mesh is regenerated at each time step allows us to make microscopic and macroscopic particles to be coincident with fluid nodes. An advantage of this procedure is that it provides an accurate definition of the particles at each time step which is useful for FSI problems.
830
831
The algorithm to compute the position of the particles using the nodal algorithm is as follows.
832
833
At each time step <math display="inline">{}^n t</math>:
834
835
<ol>
836
837
<li>Compute the converged value of the position of the fluid nodes (<math display="inline">{}^{n+1}{x}_i,~i = 1,\cdots ,N</math>) and the particles (<math display="inline">{}^{n+1}{x}_j,~j = 1,\cdots ,N_p</math>) using the algorithm of Section 9. The <math display="inline">N_p</math> particles coinciding with <math display="inline">N_p</math> fluid nodes (<math display="inline">N_p \le N</math>) will typically move to a different position than that of the corresponding fluid nodes (Figure [[#img-9|9]]).   </li>
838
<li>Regenerate the mesh discretizing the fluid domain treating the position of the <math display="inline">N_p</math> particles as fluid nodes and ignore the fluid nodes originally coinciding with the <math display="inline">N_p</math> particles at <math display="inline">{}^{n+1} t</math>.   </li>
839
<li>Interpolate the velocity of the fluid nodes at the position of the <math display="inline">N_p</math> particles surrounding the fluid nodes. </li>
840
841
</ol>
842
843
The algorithm is schematically described in Figure [[#img-9|9]].
844
845
Figures [[#img-10|10]] show an example of the application of the nodal algorithm to the study of  the motion of an individual particle  within a rectangular domain filled with water. The correct end velocity for the individual particle is obtained as shown in Figure [[#img-10|10]]c. Figures [[#img-11|11]]&#8211;[[#img-13|13]] show other examples of application of the nodal algorithm to the motion of macro-particles in water containers.
846
847
Other examples of application of this algorithm are shown in the next section.
848
849
<div id='img-9'></div>
850
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
851
|-
852
|[[Image:draft_Samper_928603926-motion-particle-node.png|600px|Nodal algorithm for tracking the motion of particles submerged in a fluid. (a) Particle i is coincident with a fluid node. (b) Update the position of the particle and the adjacent nodes. (c) Regeneration of the fluid mesh consistent with the new particle position]]
853
|- style="text-align: center; font-size: 75%;"
854
| colspan="1" | '''Figure 9:''' Nodal algorithm for tracking the motion of particles submerged in a fluid. (a) Particle <math>i</math> is coincident with a fluid node. (b) Update the position of the particle and the adjacent nodes. (c) Regeneration of the fluid mesh consistent with the new particle position
855
|}
856
857
<div id='img-10'></div>
858
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
859
|-
860
|[[Image:draft_Samper_928603926-one-particle1.png|400px|]]
861
|[[Image:draft_Samper_928603926-one-particle-contours.png|400px|]]
862
|-
863
| colspan="2"|[[Image:draft_Samper_928603926-Grafica-velocity.png|360px|Cylindrical particles falling in a water container. 2D PFEM solution using the nodal algorithm for tracking the particle motion. (a) Mesh and particle at a certain instant. (b) Contours of the vertical velocity module. (c) Evolution of the vertical velocity of the particle until a steady state solution is found <span id='citeF-6'></span><span id='citeF-15'></span>[[#cite-6|[6,15]]]]]
864
|- style="text-align: center; font-size: 75%;"
865
| colspan="2" | '''Figure 10:''' Cylindrical particles falling in a water container. 2D PFEM solution using the nodal algorithm for tracking the particle motion. (a) Mesh and particle at a certain instant. (b) Contours of the vertical velocity module. (c) Evolution of the vertical velocity of the particle until a steady state solution is found <span id='citeF-6'></span><span id='citeF-15'></span>[[#cite-6|[6,15]]]
866
|}
867
868
<div id='img-11'></div>
869
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
870
|-
871
|[[Image:draft_Samper_928603926-particles-crossing1.png|400px|]]
872
|[[Image:draft_Samper_928603926-particles-crossing2.png|400px|Motion of ascending and descending particles of different density in a fluid domain. PFEM results using the nodal algorithm for tracking the particles motion]]
873
|- style="text-align: center; font-size: 75%;"
874
| colspan="2" | '''Figure 11:''' Motion of ascending and descending particles of different density in a fluid domain. PFEM results using the nodal algorithm for tracking the particles motion
875
|}
876
877
<div id='img-12'></div>
878
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
879
|-
880
|[[Image:draft_Samper_928603926-three-particles2.png|600px|Motion of three macroscopic particles in a water sloshing problem within a tank. PFEM results obtained using the nodal algorithm for particle tracking]]
881
|- style="text-align: center; font-size: 75%;"
882
| colspan="1" | '''Figure 12:''' Motion of three macroscopic particles in a water sloshing problem within a tank. PFEM results obtained using the nodal algorithm for particle tracking
883
|}
884
885
<div id='img-13'></div>
886
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
887
|-
888
|[[Image:draft_Samper_928603926-tanque-particles-entrada_clara.png|400px|]]
889
|[[Image:draft_Samper_928603926-tanque-particles-clara_clara.png|400px|PFEM analysis of the penetration of a collection of spherical (macroscopic) particles into a water container]]
890
|- style="text-align: center; font-size: 75%;"
891
| colspan="2" | '''Figure 13:''' PFEM analysis of the penetration of a collection of spherical (macroscopic) particles into a water container
892
|}
893
894
<div id='img-14'></div>
895
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
896
|-
897
|[[Image:draft_Samper_928603926-Fig-soilmass1.png|400px|]]
898
|[[Image:draft_Samper_928603926-Fig-soilmass2.png|400px|]]
899
|-
900
|[[Image:draft_Samper_928603926-Fig-soilmass3.png|400px|]]
901
|[[Image:draft_Samper_928603926-Fig-soilmass4.png|400px|Falling of a lorry into the sea by erosion of the underlying soil mass due to the action of waves]]
902
|- style="text-align: center; font-size: 75%;"
903
| colspan="2" | '''Figure 14:''' Falling of a lorry into the sea by erosion of the underlying soil mass due to the action of waves
904
|}
905
906
<div id='img-15'></div>
907
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
908
|-
909
|[[Image:draft_Samper_928603926-caida-particulas1.png|400px|]]
910
|[[Image:draft_Samper_928603926-caida-particulas2.png|400px|Sliding of macroscopic particles over an inclined wall entering a container with water]]
911
|- style="text-align: center; font-size: 75%;"
912
| colspan="2" | '''Figure 15:''' Sliding of macroscopic particles over an inclined wall entering a container with water
913
|}
914
915
<div id='img-16'></div>
916
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
917
|-
918
|[[Image:draft_Samper_928603926-landslide-falling1.png|400px|]]
919
|[[Image:draft_Samper_928603926-landslide-falling2.png|400px|3D PFEM simulation of a landslide falling on four houses]]
920
|- style="text-align: center; font-size: 75%;"
921
| colspan="2" | '''Figure 16:''' 3D PFEM simulation of a landslide falling on four houses
922
|}
923
924
<div id='img-17'></div>
925
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
926
|-
927
|[[Image:draft_Samper_928603926-picture_1_of_4.png|400px|]]
928
|[[Image:draft_Samper_928603926-picture_2_of_4.png|400px|]]
929
|-
930
|[[Image:draft_Samper_928603926-picture_3_of_4.png|400px|]]
931
|[[Image:draft_Samper_928603926-picture_4_of_4.png|400px|3D PFEM results for the dragging of a collection of large rocks by a water stream]]
932
|- style="text-align: center; font-size: 75%;"
933
| colspan="2" | '''Figure 17:''' 3D PFEM results for the dragging of a collection of large rocks by a water stream
934
|}
935
936
==12 Examples==
937
938
We present the study of a several problems involving the transport of macroscopic and large particles in fluid flows. The problems are schematic representations of particulate flows occurring in practical problems of civil and environmental engineering and industrial problems.
939
940
Most problems studied have been solved with the PFEM using the  nodal algorithm for the transport of macroscopic particles described in the previous section. An exception are the problems in Section 12.6 dealing with the vertical  transport of spherical particles in a cylinder where  the standard immersed  approach for the transport macroparticles described in Sections 1&#8211;4 was used and the fluid equations were solved with an Eulerian flow solver implemented in the Kratos open-source software platform of CIMNE <span id='citeF-24'></span>[[#cite-24|[24]]].
941
942
===12.1 Erosion of a slope adjacent to the shore due to sea waves===
943
944
Figure [[#img-14|14]] shows a representative example of the progressive erosion of a soil mass adjacent to the shore due to sea waves and the subsequent falling into the sea of a 2D object representing the section of a lorry. The object has been modeled as a rigid solid. Note that the eroded soil particles accumulate at the sea bottom.
945
946
This example, although still quite simple and schematic, evidences the possibility of the PFEM  for modeling FSSI problems involving soil erosion, free surface waves and rigid/deformable structures.
947
948
===12.2 Landslide falling on houses===
949
950
Figure [[#img-15|15]] shows two instants of the 2D simulation with PFEM of the motion of a collection of macroscopic particles as they slide over an inclined wall and fall into a water container.
951
952
The PFEM is particularly suited for the modelling and simulation of landslides and their effect in the surrounding structures. Figure [[#img-16|16]] shows an schematic 2D simulation of a landslide falling on two adjacent constructions. The landslide material has been assumed to behave as a pure viscoplastic material modelled as a non-Newtonian viscous incompressible fluid <span id='citeF-57'></span>[[#cite-57|[57]]]. Other applications of the PFEM to the modelling of landslides can be found in <span id='citeF-8'></span><span id='citeF-49'></span>[[#cite-8|[8,49]]].
953
954
===12.3 Dragging of rocks by a water stream===
955
956
Figure [[#img-17|17]] shows the dragging of a collection of rocks modelled as large rigid particles of arbitrary shape by the action of a water stream. The particles move due to the action of the water forces on the particles computed by integrating the pressure and the viscous stresses in the elements surrounding each particle.
957
958
===12.4 Suction of cohesive material submerged in water===
959
960
Figures [[#img-18|18]] and [[#img-19|19]] show two examples of the detachment, suction and transport of particles of a cohesive material submerged on water. Figure [[#img-18|18]] shows how the particles detatch  from the cohesive soil bed and are transported within the suctioning tube (modelled as a 2D problem). The last picture shows the erosion in the soil as the mixture of water and eroded particles falls down from within the tube and hits the soil surface due to a stop in the suction pressure.
961
962
Figure [[#img-19|19]] shows a similar 3D problem. The suction in the tube erodes the surface of the soil bed in the right hand container. The mixture of water and eroded particles is transported to the adjacent containers.
963
964
<div id='img-18'></div>
965
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
966
|-
967
|[[Image:draft_Samper_928603926-suction1.png|400px|]]
968
|[[Image:draft_Samper_928603926-suction2.png|400px|]]
969
|-
970
|[[Image:draft_Samper_928603926-suction3.png|400px|]]
971
|[[Image:draft_Samper_928603926-suction4.png|400px|2D PFEM analysis of the detachment and suction of cohesive material submerged in water. The last picture shows the erosion of the bed material after the impact of the mixture of water and eroded particles falling from within  the tube]]
972
|- style="text-align: center; font-size: 75%;"
973
| colspan="2" | '''Figure 18:''' 2D PFEM analysis of the detachment and suction of cohesive material submerged in water. The last picture shows the erosion of the bed material after the impact of the mixture of water and eroded particles falling from within  the tube
974
|}
975
976
<div id='img-19'></div>
977
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
978
|-
979
|[[Image:draft_Samper_928603926-picture1_of_4.png|400px|]]
980
|[[Image:draft_Samper_928603926-picture2_of_4.png|400px|]]
981
|-
982
|[[Image:draft_Samper_928603926-picture3_of_4.png|400px|]]
983
|[[Image:draft_Samper_928603926-picture4_of_4.png|400px|3D PFEM simulation of the detachment, suction and transport of submerged cohesive material from one recipient to another]]
984
|- style="text-align: center; font-size: 75%;"
985
| colspan="2" | '''Figure 19:''' 3D PFEM simulation of the detachment, suction and transport of submerged cohesive material from one recipient to another
986
|}
987
988
===12.5 Filling of a water container with particles===
989
990
Figure [[#img-20|20]] shows a 3D example of the filling of a cylindrical container  with water containing macroscopic spherical particles. Water is allowed to exit the cylinder by a lateral hole while particles enter from two other holes at different height and fall down by gravity until they progressively fill the cylinder. The figures show different instants of the filling process.
991
992
<div id='img-20'></div>
993
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
994
|-
995
|[[Image:draft_Samper_928603926-diposit_2_of_6.png|400px|]]
996
|[[Image:draft_Samper_928603926-diposit_3_of_6.png|400px|]]
997
|-
998
|[[Image:draft_Samper_928603926-diposit_5_of_6.png|400px|]]
999
|[[Image:draft_Samper_928603926-diposit_6_of_6.png|400px|Filling of a container by injecting water containing macroscopic particles from two holes. Water is allowed to exit through a third hole at the upper right hand side of the cylinder. 3D PFEM results at four instants ]]
1000
|- style="text-align: center; font-size: 75%;"
1001
| colspan="2" | '''Figure 20:''' Filling of a container by injecting water containing macroscopic particles from two holes. Water is allowed to exit through a third hole at the upper right hand side of the cylinder. 3D PFEM results at four instants 
1002
|}
1003
1004
===12.6 Transport of spherical particles in a tube filled with water===
1005
1006
The example in Figure [[#img-21|21]] models the vertical transport of some 120.000 spherical particles to the surface of a tube filled with water and the subsequent deposition of the particles on the free water surface at the top of the tube. Particles move upwards within the tube due to a fluid velocity of 1 m/s. The average size of the particles radius is 2 mm and their density is 2300 Kg/m<math display="inline">^3</math>. Particles move vertically until they reach the top of the fluid domain and accumulate there due to the combined effect of their weight and the effect of the interaction force from the fluid. Figure [[#img-21|21]] shows two instants of the particles ascending process. The accumulation of particles in the water free surface  at the  top of the tube is clearly seen.
1007
1008
Figure [[#img-22|22]] shows the interaction of eigth jets of ascending air bubbles with 200.000 spherical solid particles that fall down within a cylinder filled with water.
1009
1010
<div id='img-21'></div>
1011
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1012
|-
1013
|[[Image:draft_Samper_928603926-ascending-particles1b.png|400px|]]
1014
|[[Image:draft_Samper_928603926-ascending-particles2.png|400px|Transport of spherical macroparticles up to the free surface of a tube filled with water. Particles move up with a prescribed velocity until they accumulate on the free surface. Results obtained with a coupled DEM-Eulerian CFD code <span id='citeF-24'></span>[[#cite-24|[24]]]]]
1015
|- style="text-align: center; font-size: 75%;"
1016
| colspan="2" | '''Figure 21:''' Transport of spherical macroparticles up to the free surface of a tube filled with water. Particles move up with a prescribed velocity until they accumulate on the free surface. Results obtained with a coupled DEM-Eulerian CFD code <span id='citeF-24'></span>[[#cite-24|[24]]]
1017
|}
1018
1019
<div id='img-22'></div>
1020
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1021
|-
1022
|[[Image:draft_Samper_928603926-interaction-six-jets1.png|400px|]]
1023
|[[Image:draft_Samper_928603926-interaction-six-jets2.png|400px|]]
1024
|-
1025
|[[Image:draft_Samper_928603926-interaction-six-jets3.png|400px|]]
1026
|[[Image:draft_Samper_928603926-interaction-six-jets4.png|400px|Interaction of eight jets of ascending air bubbles with a thick layer of 200.000 spherical particles that fall down within a cylinder filled with water. Numerical results obtained with a coupled DEM-Eulerian CFD code <span id='citeF-24'></span>[[#cite-24|[24]]]]]
1027
|- style="text-align: center; font-size: 75%;"
1028
| colspan="2" | '''Figure 22:''' Interaction of eight jets of ascending air bubbles with a thick layer of 200.000 spherical particles that fall down within a cylinder filled with water. Numerical results obtained with a coupled DEM-Eulerian CFD code <span id='citeF-24'></span>[[#cite-24|[24]]]
1029
|}
1030
1031
===12.7 Dragging of large objects and small particles in a tsunami type flow===
1032
1033
The last example is the dragging of cars, barrels and debris  (modelled as macroscopic particles) by a water stream that flows over a vertical wall. The problem represents an schematic study of a real situation corresponding to the tsunami in Fukushima (Japan) on March 2011 (Figure [[#img-23|23]]). Figures [[#img-24|24]] show two snapshots of the PFEM solution of this complex problem.
1034
1035
<div id='img-23'></div>
1036
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
1037
|-
1038
|[[Image:draft_Samper_928603926-tsunami-japon.png|450px|Dragging of cars and large and small objects in the Fukushima tsunami (Japan) ]]
1039
|- style="text-align: center; font-size: 75%;"
1040
| colspan="1" | '''Figure 23:''' Dragging of cars and large and small objects in the Fukushima tsunami (Japan) 
1041
|}
1042
1043
<div id='img-24'></div>
1044
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1045
|-
1046
|[[Image:draft_Samper_928603926-vlcsnap-2013-05-27-13h19m40s59.png|600px|]]
1047
|-
1048
|[[Image:draft_Samper_928603926-vlcsnap-2013-05-27-13h21m48s219.png|510px|Dragging of large objects and macroscopic particles in a tsunami type flow passing over a vertical wall. 3D PFEM results ]]
1049
|- style="text-align: center; font-size: 75%;"
1050
| colspan="1" | '''Figure 24:''' Dragging of large objects and macroscopic particles in a tsunami type flow passing over a vertical wall. 3D PFEM results 
1051
|}
1052
1053
==13 Concluding remarks==
1054
1055
We have presented a Lagrangian numerical technique for analysis of flows incorporating physical particles of different sizes. The numerical approach is based on the Particle Finite Element Method (PFEM) and a stabilized Lagrangian mixed velocity-pressure formulation. The examples presented in the paper evidence the possibilities of the PFEM for analysis of practical multiscale particulate flows in industrial and environmental problems.
1056
1057
==Acknowledgements==
1058
1059
This research was  supported by the Advanced Grant  project SAFECON of the European Research Council.
1060
1061
===BIBLIOGRAPHY===
1062
1063
<div id="cite-1"></div>
1064
'''[[#citeF-1|[1]]]'''  Anderson T, Jackson R. Fluid mechanical description of fluidized beds: equation of motion. Industrial & Engineering Chemical Fundamentals, 6(4):527&#8211;539
1065
1066
<div id="cite-2"></div>
1067
'''[[#citeF-2|[2]]]'''  Avci B, Wriggers P, (2012) A DEM-FEM coupling approach for the direct numerical simulation of 3D particulate flows. Journal of Applied Mechanics, 79(1), 7 pages
1068
1069
<div id="cite-3"></div>
1070
'''[[#citeF-3|[3]]]'''  T. Belytschko, W.K. Liu, B. Moran, Non linear finite element for continua and structures, 2d Edition, Wiley, 2013.
1071
1072
<div id="cite-4"></div>
1073
'''[[#citeF-4|[4]]]'''  Carbonell JM, Oñate E, Suárez B (2010) Modeling of ground excavation with the Particle Finite Element Method. Journal of Engineering Mechanics (ASCE) 136(4):455&#8211;463
1074
1075
<div id="cite-5"></div>
1076
'''[[#citeF-5|[5]]]'''  Carbonell JM, Oñate E, Suárez B (2013) Modelling of tunnelling processes and cutting tool wear with the Particle Finite Element Method (PFEM). Accepted in Comput. Mech. (2013) DOI:10.1007/s00466-013-0835-x
1077
1078
<div id="cite-6"></div>
1079
'''[[#citeF-6|[6]]]'''  Clift R, Grace JR, Weber ME (1978) Bubbles, drops and particles. Academic Press, New York
1080
1081
<div id="cite-7"></div>
1082
'''[7]'''  Coussy O (2004) Poromechanics. Wiley
1083
1084
<div id="cite-8"></div>
1085
'''[[#citeF-8|[8]]]'''  Cremonesi M, Frangi A, Perego U (2011) A Lagrangian finite element approach for the simulation of water-waves induced by landslides. Computers & Structures 89:1086&#8211;1093
1086
1087
<div id="cite-9"></div>
1088
'''[[#citeF-9|[9]]]'''   Donea J,  Huerta A  (2003) Finite element method for flow problems. J. Wiley
1089
1090
<div id="cite-10"></div>
1091
'''[[#citeF-10|[10]]]'''  Edelsbrunner H,  Mucke EP (1999) Three dimensional alpha shapes. ACM   Trans. Graphics  13:43&#8211;72
1092
1093
<div id="cite-11"></div>
1094
'''[[#citeF-11|[11]]]'''  Franci A, Oñate E, Carbonell JM (2013) On the effect of the  tangent bulk stiffness matrix in the analysis of free surface Lagrangian flows using PFEM. Research Report CIMNE PI402
1095
1096
<div id="cite-12"></div>
1097
'''[[#citeF-12|[12]]]'''  Franci A, Oñate E, Carbonell JM (2014b) Unified updated Lagrangian formulation for fluid-structure interaction problems. Research Report CIMNE PI404
1098
1099
<div id="cite-13"></div>
1100
'''[[#citeF-13|[13]]]'''  Gidaspow D (1994) Multiphase flow and fluidization. Continuum and Kinetic Theory Description, Academic Press, 467 pages
1101
1102
<div id="cite-14"></div>
1103
'''[[#citeF-14|[14]]]'''  Healy DP, Young DB (2005) Full Lagrangian method for calculating particle concentration field in dilute gas-particle flows. Proc. Roy. Soc., London A: Mathematical, Physical and Engineering Sciences, 461(2059):2197&#8211;2225
1104
1105
<div id="cite-15"></div>
1106
'''[[#citeF-15|[15]]]'''  Heider A, Levespiel O (1989) Drag coefficient and terminal velocity of spherical and nonspherical particles. Powder Technol. 58:63&#8211;70
1107
1108
<div id="cite-16"></div>
1109
'''[[#citeF-16|[16]]]''' Hilton J, Cleary P (2013) Dust modelling using a combined CFD and discrete element formulation. Int. J. Num. Meth. Fluids, 72(5):528-549
1110
1111
<div id="cite-17"></div>
1112
'''[[#citeF-17|[17]]]'''  Idelsohn SR, Oñate E, Del Pin F (2004)  The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int. Journal for Numerical Methods in Engineering,  61(7):964&#8211;989
1113
1114
<div id="cite-18"></div>
1115
'''[[#citeF-18|[18]]]'''  Idelsohn SR, Marti J, Limache A, Oñate E (2008) Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluid-structure interaction problems via the PFEM. Comput Methods Appl Mech Engrg.  197:1762&#8211;1776
1116
1117
<div id="cite-19"></div>
1118
'''[19]'''  Idelsohn SR, Mier-Torrecilla M, Oñate E  (2009) Multi-fluid flows with the Particle Finite Element Method. Comput Methods Appl Mech Engrg. 198:2750&#8211;2767
1119
1120
<div id="cite-20"></div>
1121
'''[[#citeF-20|[20]]]'''  Idelsohn SR,  Oñate E (2010) The challenge of mass conservation in the solution of free-surface flows with the fractional-step method: Problems and solutions. Int. J. Numer. Meth. Biomed. Engng. 26:1313–-1330
1122
1123
<div id="cite-21"></div>
1124
'''[[#citeF-21|[21]]]'''  Jajcevic D, Siegmann E, Radeke C, Khinast JG (2013) Large-scale CFD-DEM simulations of fluidized granular systems. Chemical Engineering Science, 98:298&#8211;310
1125
1126
<div id="cite-22"></div>
1127
'''[[#citeF-22|[22]]]'''  Jackson R (2000), The dynamics of fluidized particles. Cambridge Monographs on Mechanics, Cambridge Univ. Press
1128
1129
<div id="cite-23"></div>
1130
'''[[#citeF-23|[23]]]'''  Kafui DK, Thornton C, Adams MJ (2002) Discrete particle-continuum fluid modelling of gas-solid fluidised beds.  Chemical Engng. Science, 57(13):2395&#8211;2410
1131
1132
<div id="cite-24"></div>
1133
'''[[#citeF-24|[24]]]'''  Kratos (2014) Open source software platform for multiphysics computations. CIMNE, Barcelona, www.cimne.com/kratos
1134
1135
<div id="cite-25"></div>
1136
'''[25]'''  Larese A,  Rossi R, Oñate E, Idelsohn SR (2008)  Validation of the Particle Finite Element Method (PFEM) for simulation of free surface flows. Engineering Computations 25(4):385&#8211;425
1137
1138
<div id="cite-26"></div>
1139
'''[[#citeF-26|[26]]]'''  Liu SH, Sun DA (2002). Simulating the collapse of unsaturated soil by DEM. Int. J. Num. Anal. Meth. in Geomechanics,  26:633–646
1140
1141
<div id="cite-27"></div>
1142
'''[[#citeF-27|[27]]]'''  Oliver X, Cante JC, Weyler R, González C, Hernández J (2007) Particle finite element methods in solid mechanics problems. In: Oñate E, Owen R (Eds) Computational Plasticity. Springer, Berlin, pp 87&#8211;103
1143
1144
<div id="cite-28"></div>
1145
'''[[#citeF-28|[28]]]'''  Oliver X, Hartmann S, Cante JC, Wyler R, Hernández J (2009) A contact domain method for large deformation frictional contact problems. Part 1: theoretical basis. Comput Methods Appl Mech Eng 198:2591&#8211;2606
1146
1147
<div id="cite-29"></div>
1148
'''[[#citeF-29|[29]]]'''  Oñate E (1998) Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Engng. 151:233&#8211;267
1149
1150
<div id="cite-30"></div>
1151
'''[[#citeF-30|[30]]]'''  Oñate E (2000) A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comput Methods Appl Mech Engrg. 182(1&#8211;2):355&#8211;370
1152
1153
<div id="cite-31"></div>
1154
'''[[#citeF-31|[31]]]'''   Oñate E, García J (2001)  A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation. Comput. Meth. Appl. Mech. Engrg. 191:635&#8211;660
1155
1156
<div id="cite-32"></div>
1157
'''[[#citeF-32|[32]]]'''  Oñate E (2003) Multiscale computational analysis in mechanics using   finite calculus: an introduction. Comput. Meth. Appl. Mech. Engrg.   192(28-30):3043&#8211;3059
1158
1159
<div id="cite-33"></div>
1160
'''[[#citeF-33|[33]]]'''   Oñate E (2004) Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng. 60(1):255&#8211;281
1161
1162
<div id="cite-34"></div>
1163
'''[[#citeF-34|[34]]]'''   Oñate E, Rojek J,  (2004b) Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems. Comput. Meth. Appl. Mech. Engrg. 193:3087&#8211;3128
1164
1165
<div id="cite-35"></div>
1166
'''[[#citeF-35|[35]]]'''  Oñate E, Idelsohn SR, Del Pin F, Aubry R (2004c) The particle   finite element method. An overview. Int. J. Comput. Methods    1(2):267&#8211;307
1167
1168
<div id="cite-36"></div>
1169
'''[[#citeF-36|[36]]]'''  Oñate E, Celigueta MA, Idelsohn SR (2006a) Modeling bed erosion in   free surface flows by the Particle Finite Element Method, Acta   Geotechnia 1(4):237&#8211;252
1170
1171
<div id="cite-37"></div>
1172
'''[[#citeF-37|[37]]]'''  Oñate E, Valls A, García J  (2006b) FIC/FEM formulation   with matrix stabilizing terms for incompressible flows at low and high   Reynold's numbers. Computational Mechanics 38 (4-5):440&#8211;455
1173
1174
<div id="cite-38"></div>
1175
'''[[#citeF-38|[38]]]'''  Oñate E, García J,  Idelsohn SR, Del Pin F  (2006c)  FIC   formulations for finite element analysis of incompressible flows. Eulerian,   ALE and Lagrangian approaches.  Comput. Meth. Appl. Mech. Engng.   195(23-24):3001&#8211;3037
1176
1177
<div id="cite-39"></div>
1178
'''[[#citeF-39|[39]]]'''  Oñate E, Valls A, García J  (2007) Computation of turbulent flows using a finite calculus-finite element formulation. Int. J. Numer. Meth. Engng. 54:609&#8211;637
1179
1180
<div id="cite-40"></div>
1181
'''[[#citeF-40|[40]]]'''  Oñate E, Idelsohn SR,  Celigueta MA, Rossi R (2008) Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. Comput. Meth. Appl. Mech. Engng. 197(19-20):1777&#8211;1800
1182
1183
<div id="cite-41"></div>
1184
'''[[#citeF-41|[41]]]'''  Oñate E (2009),  Structural analysis with the finite element method. Linear statics. Volume 1. Basis and Solids.  CIMNE-Springer
1185
1186
<div id="cite-42"></div>
1187
'''[42]'''  Oñate E, Rossi R,  Idelsohn SR, Butler K (2010) Melting and spread of polymers in fire with the particle finite element method. Int. J. Numerical Methods in Engineering, 81(8):1046&#8211;1072
1188
1189
<div id="cite-43"></div>
1190
'''[[#citeF-43|[43]]]'''  Oñate E, Celigueta MA, Idelsohn SR,  Salazar F, Suárez B (2011) Possibilities of the particle finite element method for fluid-soil-structure interaction problems. Computational Mechanics, 48(3):307&#8211;318.
1191
1192
<div id="cite-44"></div>
1193
'''[[#citeF-44|[44]]]'''  Oñate E, Carbonell JM (2013)  Updated Lagrangian finite element formulation for quasi and fully incompressible fluids. Research Report PI393, CIMNE. Submitted to Comput. Mechanics
1194
1195
<div id="cite-45"></div>
1196
'''[[#citeF-45|[45]]]'''  Oñate E, Franci A, Carbonell JM (2014) Lagrangian formulation  for finite element analysis of quasi-incompressible fluids with reduced mass losses. Accepted for publication in Int. J. Num. Meth. in Fluids, 74(10):699–-731
1197
1198
<div id="cite-46"></div>
1199
'''[[#citeF-46|[46]]]'''  Oñate E, Franci A, Carbonell JM (2014b) A Particle Finite Element Method (PFEM) for analysis of industrial forming processes. Accepted for publication in Comput. Mechanics
1200
1201
<div id="cite-47"></div>
1202
'''[[#citeF-47|[47]]]'''  Patankar NA, Joseph DD (2001) Lagrangian numerical simulation of particulate flows. Int. J. Multiphase Flow, 27:1685&#8211;1706
1203
1204
<div id="cite-48"></div>
1205
'''[48]'''  Ryzhakov P, Oñate E, Rossi R, Idelsohn SR (2012) Improving mass conservation in simulation of incompressible flows. Int. J. Numer. Meth. Engng.  90(12):1435&#8211;1451
1206
1207
<div id="cite-49"></div>
1208
'''[[#citeF-49|[49]]]'''  Salazar F, Oñate E, Morán R (2012) Modelación numérica del deslizamiento  de ladeu en embalses mediante el método de partículos y elementos finitos (PFEM). Rev. Int. Meth. Num. Cal. Dis. Ing., 28(2):112&#8211;123
1209
1210
<div id="cite-50"></div>
1211
'''[[#citeF-50|[50]]]'''  Shamy UE, Zeghal M (2005)  Coupled continuum&#8211;discrete model for saturated granular soils. J. Engineering Mechanics (ASCE), 131(4):413–-426
1212
1213
<div id="cite-51"></div>
1214
'''[[#citeF-51|[51]]]'''  Sommerfeld M, van Wachen B, Oliemans R (Eds) (2008) Best practice guideliens for computational fluid dynamics of dispersed multiphase flows. European Research Community on Flow, Turbulence and Combustion (ERCOFTAC), Imperial College London.
1215
1216
<div id="cite-52"></div>
1217
'''[52]'''  Tang B, Li JF, Wang TS (2009) Some improvements on free surface simulation by the particle finite element method. Int. J. Num. Methods in Fluids, 60(9):1032–-1054
1218
1219
<div id="cite-53"></div>
1220
'''[[#citeF-53|[53]]]'''  van Wachen B, Oliveira ES (2010) An immersed boundary method for interacting particles. ERCOFTAC Bulletin 82, 17&#8211;22 March 2010
1221
1222
<div id="cite-54"></div>
1223
'''[[#citeF-54|[54]]]'''  Wang X, Zhang LT, Liu WK (2009) On computational issues of the immersed finite element method. J. Comp. Physics, 228:2535&#8211;2551
1224
1225
<div id="cite-55"></div>
1226
'''[[#citeF-55|[55]]]'''  Li X, Chu X, Sheng DC (2007) A saturated discrete particle model and characteristic-based SPH method in granular materials. Int. J. Numer. Meth. Engng., 72:858&#8211;882
1227
1228
<div id="cite-56"></div>
1229
'''[[#citeF-56|[56]]]''' Zhang LT, Gerstenberg A, Wang X, Liu WK (2004) Immersed finite element method. Comput. Meth. Appl. Mech. Engrg., 193(21-22):2051&#8211;2067
1230
1231
<div id="cite-57"></div>
1232
'''[[#citeF-57|[57]]]'''  Zienkiewicz OC, Jain PC, Oñate E (1978) Flow of solids during forming and extrusion: some aspects of numerical solutions. ''Int. J. Solids Struct.'', 14:15&#8211;38
1233
1234
<div id="cite-58"></div>
1235
'''[[#citeF-58|[58]]]'''  Zienkiewicz OC,  Taylor RL, Zhu JZ (2005) The finite element method. The basis.  6th Ed., Elsevier
1236
1237
<div id="cite-59"></div>
1238
'''[[#citeF-59|[59]]]'''  Zienkiewicz OC,  Taylor RL (2005) The finite element method for   solid and structural mechanics.  6th Ed., Elsevier
1239
1240
<div id="cite-60"></div>
1241
'''[[#citeF-60|[60]]]'''  Zienkiewicz OC,  Taylor RL, Nithiarasu P (2005)  The finite element   method for fluid dynamics.  6th Ed.,  Elsevier
1242
1243
<div id="cite-61"></div>
1244
'''[[#citeF-61|[61]]]'''  Zohdi T (2007) An introduction to modelling and simulation of particulate flows. SIAM, Computational Science and Engineering
1245
1246
<div id="cite-62"></div>
1247
'''[[#citeF-62|[62]]]'''  Zohdi T, Wriggers P (2007) Computation of strongly coupled multifield interaction in particle-fluid systems. Comput. Meth. Appl. Mech. Engng., 196:3927&#8211;3950
1248

Return to Onate et al 2014d.

Back to Top

Document information

Published on 01/01/2014

DOI: 10.1007/s40571-014-0012-9
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 48
Views 40
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?