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
Published in''Numerical Simulations of Coupled problems in Engineering'', S.R. Idelsohn (Ed.). Computational Methods in Applied Sciences 33, pp. 129-156, Springer 2014<br />
2
DOI: 10.1007/978-3-319-06136-8
3
4
5
6
==Abstract==
7
8
We present a Lagrangian formulation for coupled thermal analysis of quasi and fully incompressible flows and fluid-structure interaction (FSI) problems that has excellent mass preservation features. The success of the formulation lays on  a  residual-based stabilized expression of the mass balance equation obtained using the Finite Calculus (FIC) method. The  governing equations are discretized with the FEM using simplicial elements with equal linear interpolation for the velocities, the pressure and the temperature. The merits of the formulation in terms of reduced mass loss and overall accuracy are verified in the solution of 2D and 3D adiabatic and thermally-coupled quasi-incompressible free-surface flow problems using the Particle Finite Element Method (PFEM). Examples include the sloshing of water in a tank and the falling  of a water sphere and a cylinder into a  tank containing water.
9
10
'''keywords''' Particle Finite Element Method, Coupled Thermal Analysis, Quasi and Fully Incompressible
11
12
==1 INTRODUCTION==
13
14
The analysis of thermally coupled flows and their interaction with structures is relevant in many fields of engineering. In this work we present a Lagrangian numerical technique for solving this kind of problems for quasi and fully incompressible fluids using the Particle Finite Element Method (PFEM, www.cimne.com/pfem).
15
16
The PFEM treats the mesh nodes in the analysis domain as particles which can freely move and even separate from the  domain representing, for instance, the effect of water drops or cutting particles in drilling problems. A mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM. Examples of application of PFEM to problems in fluid and solid mechanics including fluid-structure interaction (FSI) situations can be found in <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-42'></span><span id='citeF-43'></span>[[#cite-4|[4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,28,29,35,36,39,40,41,42,43]]]. Early attempts of the PFEM for solving thermally coupled flows were reported in <span id='citeF-1'></span><span id='citeF-2'></span>[[#cite-1|[1,2]]].
17
18
In Lagrangian analysis procedures (such as  PFEM) the motion of the fluid particles is tracked during the transient solution. Hence, the convective terms vanish in the momentum and heat transfer equations and no numerical stabilization is needed for treating those terms. Two other sources of mass loss, however, remain in the numerical solution of Lagrangian flows, i.e. that due to the treatment of the incompressibility constraint by a stabilized numerical method, and that induced by the inaccuracies in tracking the flow particles and, in particular, the free surface.
19
20
In this work the PFEM equations for analysis of thermally coupled flows and FSI problems are derived using the stabilized formulation based in the Finite Calculus (FIC) method proposed by Oñate et al. <span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-30'></span><span id='citeF-31'></span><span id='citeF-32'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-20|[20,21,22,23,24,25,26,27,30,31,32,37,38,39]]] that has excellent mass preservation features.
21
22
The lay-out of the paper is the following. In the next section we present the basic equations for conservation of linear momentum, mass and heat transfer for a quasi-incompressible fluid in a Lagrangian framework. A full incompressible fluid can be considered as a particular limit case of the former. 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,  the pressure and the temperature is presented and the relevant matrices and vectors of the discretized problem are given.  Details of the implicit solution of the Lagrangian FEM equations in time using a Newton iterative scheme are presented. The relevance of the bulk stiffness terms in the tangent matrix  for enhancing the convergence and accuracy of the iterative solution scheme is discussed. The basic steps of the PFEM for solving coupled free-surface FSI problems are described.
23
24
The efficiency and accuracy of the PFEM  technique are verified by solving a set of adiabatic and thermally coupled quasi-incompressible free surface flow problems in two (2D) and three (3D) dimensions with the PFEM. The adiabatic problems are the sloshing of water in a tank and the penetration of a water sphere into a cylindrical tank containing water. The thermally coupled problems considered are the extended 2D version of the adiabatic cases. The excellent performance of the numerical method proposed in terms of mass conservation and general accuracy is highlighted.
25
26
==2 GOVERNING EQUATIONS==
27
28
We write the governing equations for a quasi-incompressible Newtonian flow problem in the Lagrangian description as follows <span id='citeF-3'></span><span id='citeF-46'></span>[[#cite-3|[3,46]]].
29
30
===Momentum equations===
31
32
<span id="eq-1"></span>
33
{| class="formulaSCP" style="width: 100%; text-align: left;" 
34
|-
35
| 
36
{| style="text-align: left; margin:auto;" 
37
|-
38
| style="text-align: center;" | <math>\rho \frac{Dv_i}{Dt}-{\partial \sigma _{ij} \over \partial x_j}-b_i=0\quad , \quad i,j=1,n_s \quad  \hbox{in } \Omega  </math>
39
|}
40
| style="width: 5px;text-align: right;" | (1)
41
|}
42
43
In Eq.([[#eq-1|1]]), <math display="inline">\Omega </math> is the analysis domain with boundary <math display="inline">\Gamma </math>, <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, <math display="inline">\rho </math> is the density, <math display="inline">n_s</math> is the number of space dimensions (i.e. <math display="inline">n_s=3</math> for 3D problems) and <math display="inline">\sigma _{ij}</math> are the Cauchy stresses that are split in the deviatoric (<math display="inline">s_{ij}</math>) and pressure (<math display="inline">p</math>) components as
44
45
<span id="eq-2"></span>
46
{| class="formulaSCP" style="width: 100%; text-align: left;" 
47
|-
48
| 
49
{| style="text-align: left; margin:auto;" 
50
|-
51
| style="text-align: center;" | <math>\sigma _{ij} = s_{ij}+ p \delta _{ij} </math>
52
|}
53
| style="width: 5px;text-align: right;" | (2)
54
|}
55
56
where <math display="inline">\delta _{ij}</math> is the Kronecker delta. Note that the pressure is assumed to be positive for a tension state. Summation of terms with repeated indices  is assumed in Eq.([[#eq-1|1]]) and in  the following, unless otherwise specified.
57
58
The relationship between the deviatoric stresses <math display="inline">s_{ij}</math> and the strain rates <math display="inline">\varepsilon _{ij}</math> has the standard form for a Newtonian fluid,
59
60
<span id="eq-3"></span>
61
{| class="formulaSCP" style="width: 100%; text-align: left;" 
62
|-
63
| 
64
{| style="text-align: left; margin:auto;" 
65
|-
66
| style="text-align: center;" | <math>s_{ij} = 2\mu \left(\varepsilon _{ij} - {1 \over 3}\varepsilon _v \delta _{ij}\right)\quad \hbox{with}\quad  \varepsilon _{ij} = {1 \over 2}\left({\partial v_i \over \partial x_j} +  {\partial v_j \over \partial x_i}  \right) </math>
67
|}
68
| style="width: 5px;text-align: right;" | (3)
69
|}
70
71
where <math display="inline">\mu </math> is the viscosity  and <math display="inline">\varepsilon _v</math> is the volumetric strain rate defined as <math display="inline">\varepsilon _v= \varepsilon _{ii} = {\partial v_i \over \partial x_i} </math>.
72
73
===Mass balance equation===
74
75
The standard mass balance equation for a quasi-incompressible fluid can be written as <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]]
76
77
<span id="eq-4"></span>
78
{| class="formulaSCP" style="width: 100%; text-align: left;" 
79
|-
80
| 
81
{| style="text-align: left; margin:auto;" 
82
|-
83
| style="text-align: center;" | <math>r_v =0 </math>
84
|}
85
| style="width: 5px;text-align: right;" | (4a)
86
|}
87
88
with
89
90
<span id="eq-4b"></span>
91
{| class="formulaSCP" style="width: 100%; text-align: left;" 
92
|-
93
| 
94
{| style="text-align: left; margin:auto;" 
95
|-
96
| style="text-align: center;" | <math>r_v:=-\frac{1}{c^2}\frac{Dp}{Dt}+\rho \varepsilon _v </math>
97
|}
98
| style="width: 5px;text-align: right;" | (4b)
99
|}
100
101
In Eq.([[#eq-5|4b]]) <math display="inline">c</math> is the speed of sound in the fluid. For a fully incompressible fluid <math display="inline">c=\infty </math> and Eq.([[#eq-4|4a]]) simplifies to the standard form, <math display="inline">\varepsilon _v =0</math>. In our work we will retain the quasi-incompressible form of <math display="inline">r_v</math> of Eq.([[#eq-5|4b]]) for convenience.
102
103
===Thermal balance===
104
105
<span id="eq-5"></span>
106
{| class="formulaSCP" style="width: 100%; text-align: left;" 
107
|-
108
| 
109
{| style="text-align: left; margin:auto;" 
110
|-
111
| style="text-align: center;" | <math>\rho c \frac{DT}{Dt} - {\partial  \over \partial x_i} \left(k {\partial T \over \partial x_i}\right)+ Q =0  \quad i=1,n_s  \quad \hbox{in }\Omega   </math>
112
|}
113
| style="width: 5px;text-align: right;" | (5)
114
|}
115
116
where <math display="inline">T</math> is the temperature, <math display="inline">c</math> is the thermal capacity, <math display="inline">k</math> is the heat conductivity and <math display="inline">Q</math> is the heat source.
117
118
===Boundary conditions===
119
120
===''Mechanical problem''===
121
122
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
123
124
<span id="eq-6"></span>
125
{| class="formulaSCP" style="width: 100%; text-align: left;" 
126
|-
127
| 
128
{| style="text-align: left; margin:auto;" 
129
|-
130
| style="text-align: center;" | <math>v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v </math>
131
|}
132
| style="width: 5px;text-align: right;" | (6)
133
|}
134
135
<span id="eq-7"></span>
136
{| class="formulaSCP" style="width: 100%; text-align: left;" 
137
|-
138
| 
139
{| style="text-align: left; margin:auto;" 
140
|-
141
| style="text-align: center;" | <math>\sigma _{ij}n_j -t_i^p =0 \quad \hbox{on }\Gamma _t \quad i,j=1,n_s  </math>
142
|}
143
| style="width: 5px;text-align: right;" | (7)
144
|}
145
146
where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math> are the prescribed velocities and prescribed tractions on <math display="inline">\Gamma _v</math> and <math display="inline">\Gamma _t</math>, 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-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].
147
148
===''Thermal problem''===
149
150
<span id="eq-8"></span>
151
<span id="eq-9"></span>
152
{| class="formulaSCP" style="width: 100%; text-align: left;" 
153
|-
154
| 
155
{| style="text-align: left; margin:auto;" 
156
|-
157
| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on }\Gamma _\phi </math>
158
| style="width: 5px;text-align: right;" | (8)
159
|-
160
| style="text-align: center;" | <math> k {\partial \phi  \over \partial n} + q_n^p =0 \quad \hbox{on }\Gamma _q  </math>
161
| style="width: 5px;text-align: right;" | (9)
162
|}
163
|}
164
165
where <math display="inline">\phi ^p</math> and <math display="inline">q_n^p</math> are the prescribed temperature and the prescribed normal heat flux at the boundaries <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math>, respectively and <math display="inline">n</math> is the direction normal to the boundary.
166
167
'''Remark 1.''' The term <math display="inline">\frac{Dv_i}{Dt}</math> in Eq.([[#eq-1|1]]) is the ''material derivative'' of the <math display="inline">i</math>th velocity component <math display="inline">v_i</math>. This term is typically computed in a Lagrangian framework as
168
169
<span id="eq-10"></span>
170
{| class="formulaSCP" style="width: 100%; text-align: left;" 
171
|-
172
| 
173
{| style="text-align: left; margin:auto;" 
174
|-
175
| style="text-align: center;" | <math>\frac{Dv_i}{Dt} = \frac{{}^{n+1} v_i - {}^n v_i }{\Delta t}     </math>
176
|}
177
| style="width: 5px;text-align: right;" | (10)
178
|}
179
180
with
181
182
<span id="eq-11"></span>
183
{| class="formulaSCP" style="width: 100%; text-align: left;" 
184
|-
185
| 
186
{| style="text-align: left; margin:auto;" 
187
|-
188
| style="text-align: center;" | <math>{}^{n+1} v_i:=v_i ({}^{n+1} {\bf x},{}^{n+1} t)\quad ,\quad {}^n v_i:=v_i ({}^n{\bf x},{}^n t)       </math>
189
|}
190
| style="width: 5px;text-align: right;" | (11)
191
|}
192
193
where <math display="inline">{}^n v_i</math> is the velocity of the material point that has the position <math display="inline">{}^n{\bf x}</math> at time <math display="inline">t={}^nt</math>, where <math display="inline">{\bf x}</math> is the coordinates vector in a fixed Cartesian  system <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].
194
195
196
==3 STABILIZED MASS BALANCE EQUATION==
197
198
In this work we will use  the ''second order FIC form'' of the mass balance equation in space for a quasi-incompressible fluid <span id='citeF-37'></span><span id='citeF-38'></span>[[#cite-37|[37,38]]], as well as the first order FIC form of the mass balance equation in time. These forms have the following expressions:
199
200
===''Second order FIC mass balance equation in space''===
201
202
<span id="eq-12a"></span>
203
{| class="formulaSCP" style="width: 100%; text-align: left;" 
204
|-
205
| 
206
{| style="text-align: left; margin:auto;" 
207
|-
208
| style="text-align: center;" | <math>r_v + \frac{h_i^2}{12} \frac{\partial ^2 r_v}{\partial x^2_i}=0\qquad \hbox{in }\Omega \qquad i=1,n_s </math>
209
|}
210
| style="width: 5px;text-align: right;" | (12a)
211
|}
212
213
===''First order FIC mass balance equation in time''===
214
215
<span id="eq-12b"></span>
216
{| class="formulaSCP" style="width: 100%; text-align: left;" 
217
|-
218
| 
219
{| style="text-align: left; margin:auto;" 
220
|-
221
| style="text-align: center;" | <math>r_v + \frac{\delta }{2} \frac{D r_v}{D t}=0 \qquad \hbox{in }\Omega  </math>
222
|}
223
| style="width: 5px;text-align: right;" | (12b)
224
|}
225
226
Eq.([[#eq-13|12a]]) is obtained by expressing the balance of mass in a rectangular domain of finite size with dimensions <math display="inline">h_1\times h_2</math> (for 2D problems), where <math display="inline">h_i</math> are arbitrary distances, and retaining up to third order terms in the Taylor series expansions used for expressing the change of mass within the balance domain.
227
228
Eq.([[#eq-14|12b]]), on the other hand, is obtained by expressing the balance of mass in a space-time domain of infinitesimal length in space and finite dimension <math display="inline">\delta </math> in time <span id='citeF-20'></span>[[#cite-20|[20]]]. The derivation of Eqs.([[#''First order FIC mass balance equation in time''|12]]) for a 1D problem are shown in <span id='citeF-41'></span>[[#cite-41|[41]]].
229
230
The FIC terms in Eqs.([[#''First order FIC mass balance equation in time''|12]]) play the role of space and time stabilization terms respectively. In the discretized problem, the space dimensions <math display="inline">h_i</math> and the time dimension <math display="inline">\delta </math> are related to characteristic element dimensions and the time step increment, respectively as it will be explained later. Note that for <math display="inline">h_i =0</math> and <math display="inline">\delta =0</math> the standard infinitesimal form of the mass balance equation, <math display="inline">r_v =0</math>, is obtained.
231
232
After some transformations the stabilized mass balance equation ([[#eq-13|12a]]) is written as <span id='citeF-41'></span>[[#cite-41|[41]]]
233
234
<span id="eq-13"></span>
235
{| class="formulaSCP" style="width: 100%; text-align: left;" 
236
|-
237
| 
238
{| style="text-align: left; margin:auto;" 
239
|-
240
| style="text-align: center;" | <math>-\frac{1}{\kappa }\frac{Dp}{Dt}+\varepsilon _v  - \frac{\tau }{c^2} \frac{D^2p}{Dt^2} +\tau {\partial \hat r_{m_i} \over \partial x_i} =0  </math>
241
|}
242
| style="width: 5px;text-align: right;" | (13)
243
|}
244
245
where <math display="inline">\tau </math> is a ''stabilization parameter'' given by <span id='citeF-41'></span>[[#cite-41|[41]]]
246
247
<span id="eq-14"></span>
248
{| class="formulaSCP" style="width: 100%; text-align: left;" 
249
|-
250
| 
251
{| style="text-align: left; margin:auto;" 
252
|-
253
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{h^2}+ \frac{2\rho }{\delta } \right)^{-1}  </math>
254
|}
255
| style="width: 5px;text-align: right;" | (14)
256
|}
257
258
and <math display="inline">\hat r_{m_i}</math> is a ''static momentum term'' defined as
259
260
<span id="eq-15"></span>
261
{| class="formulaSCP" style="width: 100%; text-align: left;" 
262
|-
263
| 
264
{| style="text-align: left; margin:auto;" 
265
|-
266
| style="text-align: center;" | <math>\hat r_{m_i}= {\partial  \over \partial x_j} (2\mu \varepsilon _{ij}) + \frac{\partial p}{\partial x_i} + b_i   </math>
267
|}
268
| style="width: 5px;text-align: right;" | (15)
269
|}
270
271
Eq.([[#eq-13|13]]) is used as the starting point for deriving the stabilized FEM formulation as explained in the following sections.
272
273
==4 VARIATIONAL EQUATIONS==
274
275
===4.1 Momentum equations===
276
277
Multiplying Eq.(1) by arbitrary test functions <math display="inline">w_i</math> with dimensions of velocity and integrating over the analysis domain <math display="inline">\Omega </math> gives the weighted residual form of the momentum equations as <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]]
278
279
<span id="eq-16"></span>
280
{| class="formulaSCP" style="width: 100%; text-align: left;" 
281
|-
282
| 
283
{| style="text-align: left; margin:auto;" 
284
|-
285
| style="text-align: center;" | <math>\int _\Omega w_i \left(\rho \frac{Dv_i}{Dt}- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0 </math>
286
|}
287
| style="width: 5px;text-align: right;" | (16)
288
|}
289
290
Integrating by parts  the term involving <math display="inline">\sigma _{ij}</math> and using the Neumann boundary conditions ([[#eq-7|7]]) yields the weak variational form of the momentum equations as
291
292
<span id="eq-17"></span>
293
{| class="formulaSCP" style="width: 100%; text-align: left;" 
294
|-
295
| 
296
{| style="text-align: left; margin:auto;" 
297
|-
298
| style="text-align: center;" | <math>\int _\Omega w_i \rho \frac{Dv_i}{Dt} d\Omega + \int _\Omega \delta \varepsilon _{ij} \sigma _{ij} d\Omega - \int _\Omega w_i b_i d\Omega - \int _{\Gamma _t} w_i t_i^p d\Gamma =0 </math>
299
|}
300
| style="width: 5px;text-align: right;" | (17)
301
|}
302
303
where <math display="inline">\delta \varepsilon _{ij}= {\partial w_i \over \partial x_j}+{\partial w_j \over \partial x_i}</math> is an arbitrary (virtual) strain rate field. Eq.([[#eq-17|17]]) is the standard form of the ''principle of virtual power''  <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].
304
305
Substituting the expression of the stresses from Eq.([[#eq-2|2]]) into ([[#eq-17|17]]) gives
306
307
<span id="eq-18"></span>
308
{| class="formulaSCP" style="width: 100%; text-align: left;" 
309
|-
310
| 
311
{| style="text-align: left; margin:auto;" 
312
|-
313
| style="text-align: center;" | <math>\int _\Omega \! w_i \rho \frac{Dv_i}{Dt} d\Omega + \!\int _\Omega \left[\delta \varepsilon _{ij} 2\mu \left(\!\varepsilon _{ij} - \frac{1}{3}\varepsilon _{v} \delta _{ij}\!\right)+ \delta \varepsilon _{v}p\right]d\Omega - \!\int _\Omega w_i b_i d\Omega - \!\int _{\Gamma _t} w_i t_i^p d\Gamma =0  </math>
314
|}
315
| style="width: 5px;text-align: right;" | (18)
316
|}
317
318
Eq.([[#eq-18|18]]) can be written in matrix form as
319
320
<span id="eq-19"></span>
321
{| class="formulaSCP" style="width: 100%; text-align: left;" 
322
|-
323
| 
324
{| style="text-align: left; margin:auto;" 
325
|-
326
| style="text-align: center;" | <math>\displaystyle  \int _\Omega {\bf w}^T \rho \frac{D{\bf v}}{Dt} d\Omega + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf D} {\boldsymbol \varepsilon } d\Omega  + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf m} p d\Omega - \int _\Omega {\bf w}^T {\bf b}d\Omega - \!\int _{\Gamma _t} {\bf w}^T {\bf t}^p d\Gamma =0  </math>
327
|}
328
| style="width: 5px;text-align: right;" | (19)
329
|}
330
331
In Eq.([[#eq-19|19]]) <math display="inline">{\bf w},{\bf v},{\boldsymbol \varepsilon }</math> and <math display="inline">\delta {\boldsymbol \varepsilon }</math> are vectors containing the test functions, the velocities, the strain rates and the virtual strain rates respectively; <math display="inline">{\bf b}</math> and <math display="inline">{\bf t}^p</math> are body force and surface traction vectors, respectively; <math display="inline">{\bf D}</math> is the viscous constitutive matrix and <math display="inline">{\bf m}</math> is an auxiliary vector. These vectors are defined as (for 3D problems)
332
333
<span id="eq-20"></span>
334
{| class="formulaSCP" style="width: 100%; text-align: left;" 
335
|-
336
| 
337
{| style="text-align: left; margin:auto;" 
338
|-
339
| style="text-align: center;" | <math>\begin{array}{l}\textbf{w} =[w_1,w_2,w_3]^T\quad ,\quad \textbf{v} =[v_1,v_2,v_3]^T \quad ,\quad \textbf{b} =[b_1,b_2,b_3]^T \quad ,\quad  \textbf{t}^p =[t_1^p,t_2^p,t_3^p]^T\\[.5cm] {\boldsymbol \varepsilon }= [\varepsilon _{11},\varepsilon _{22}, \varepsilon _{33},\varepsilon _{12},\varepsilon _{13}, \varepsilon _{23}]^T \quad ,\quad  \delta {\boldsymbol \varepsilon }= [\delta \varepsilon _{11},\delta \varepsilon _{22}, \delta \varepsilon _{33},\delta \varepsilon _{12},\delta \varepsilon _{13}, \delta \varepsilon _{23}]^T\\[.5cm]  {\bf D}=\mu \left[                      \begin{array}{cccccc}4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\                         &  4/3 & -2/3 &0  &0  & 0 \\                         &  & 4/3 & 0 & 0 & 0 \\                         &  &  & 2 & 0 & 0 \\                        {Sym.} &  &  &  & 2 & 0 \\                         &  &  &  &  & 2 \\                      \end{array} \right]\qquad ,\qquad {\bf m}=[1,1,1,0,0,0]^T \end{array} </math>
340
|}
341
| style="width: 5px;text-align: right;" | (20)
342
|}
343
344
===4.2 Mass balance equations===
345
346
We multiply Eq.([[#eq-13|13]]) by arbitrary (continuous) test functions <math display="inline">q</math> (with dimensions of pressure) defined over the analysis domain <math display="inline">\Omega </math>. Integrating over <math display="inline">\Omega </math> gives
347
348
<span id="eq-21"></span>
349
{| class="formulaSCP" style="width: 100%; text-align: left;" 
350
|-
351
| 
352
{| style="text-align: left; margin:auto;" 
353
|-
354
| style="text-align: center;" | <math>\int _\Omega -\frac{q}{\kappa } \frac{Dp}{Dt}d\Omega -\int _\Omega q \frac{\tau }{c^2} \frac{D^2p}{Dt^2}d\Omega + \int _\Omega q \varepsilon _v d\Omega  + \int _\Omega q \tau {\partial \hat r_{m_i} \over \partial x_i} d\Omega =0  </math>
355
|}
356
| style="width: 5px;text-align: right;" | (21)
357
|}
358
359
Integrating by parts the last integral in Eq.([[#eq-21|21]]) and using ([[#eq-15|15]]) gives after some transformations <span id='citeF-39'></span><span id='citeF-40'></span>[[#cite-39|[39,40]]]
360
361
<span id="eq-22"></span>
362
{| class="formulaSCP" style="width: 100%; text-align: left;" 
363
|-
364
| 
365
{| style="text-align: left; margin:auto;" 
366
|-
367
| style="text-align: center;" | <math>\begin{array}{r}\displaystyle \int _\Omega \frac{q}{\kappa } \frac{Dp}{Dt}d\Omega + \int _\Omega q \frac{\tau }{c^2} \frac{D^2p}{Dt^2}d\Omega - \int _\Omega q \varepsilon _v d\Omega + \int _\Omega \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)d\Omega \\ \displaystyle - \int _{\Gamma _t} q \tau \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} \left(2\mu {\partial v_n \over \partial n} + p -t_n\right)\right]d\Gamma =0 \end{array}   </math>
368
|}
369
| style="width: 5px;text-align: right;" | (22)
370
|}
371
372
Expression ([[#eq-24|24]]) holds for 2D and 3D problems. The terms involving the first and second material time derivative of the pressure and the boundary term in Eq.([[#eq-24|24]]) are important for preserving the conservation of mass in  free-surface flow problems <span id='citeF-10'></span><span id='citeF-41'></span>[[#cite-10|[10,41]]].
373
374
===4.3 Thermal balance equation===
375
376
Application of the standard weighted residual method to the heat balance equations ([[#eq-5|5]]) and ([[#eq-9|9]]) leads, after standard operations, to <span id='citeF-7'></span><span id='citeF-44'></span>[[#cite-7|[7,44]]]
377
378
<span id="eq-23"></span>
379
{| class="formulaSCP" style="width: 100%; text-align: left;" 
380
|-
381
| 
382
{| style="text-align: left; margin:auto;" 
383
|-
384
| style="text-align: center;" | <math>\int _\Omega \hat w \rho c {\partial T \over \partial t} d\Omega + \int _\Omega {\partial \hat w  \over \partial x_i}  k {\partial T \over \partial x_i}d\Omega - \int _\Omega \hat w Q d\Omega + \int _{\Gamma _q} \hat w q_n^p d\Gamma =0 </math>
385
|}
386
| style="width: 5px;text-align: right;" | (23)
387
|}
388
389
where <math display="inline">\hat w</math> are the space weighting functions for the temperature.
390
391
==5 FEM DISCRETIZATION==
392
393
We discretize  the analysis domain 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> (<math display="inline">i=1,n</math>) of element <math display="inline">e</math> <span id='citeF-34'></span><span id='citeF-44'></span>[[#cite-34|[34,44]]]. The velocity components, the  pressure and the temperature 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-34'></span><span id='citeF-44'></span><span id='citeF-46'></span>[[#cite-34|[34,44,46]]]. In matrix form
394
395
<span id="eq-24"></span>
396
{| class="formulaSCP" style="width: 100%; text-align: left;" 
397
|-
398
| 
399
{| style="text-align: left; margin:auto;" 
400
|-
401
| style="text-align: center;" | <math>{\bf v} = {\bf N}_v\bar {\bf v}  ~,~p = {\bf N}_p\bar {\bf p}\quad  , \quad T = {\bf N}_T \bar {\bf T}  </math>
402
|}
403
| style="width: 5px;text-align: right;" | (24)
404
|}
405
406
where
407
408
<span id="eq-25"></span>
409
{| class="formulaSCP" style="width: 100%; text-align: left;" 
410
|-
411
| 
412
{| style="text-align: left; margin:auto;" 
413
|-
414
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \bar {\bf v} = \left\{\begin{matrix}\bar{\bf v}^1\\\bar{\bf v}^2\\\vdots \\ \bar{\bf v}^N\end{matrix}  \right\}\quad \hbox{with}~ \bar{\bf v}^i = \left\{\begin{matrix}\bar{v}^i_1\\\bar{v}^i_2\\\bar{v}^i_3\end{matrix}  \right\}~, ~ \bar {\bf p} =\left\{\begin{matrix}\bar{p}^1\\\bar{p}^2\\\vdots \\ \bar{p}^{N}\end{matrix}  \right\}~, ~ {\bf T}= \left\{\begin{matrix}\bar{T}_1\\\bar{T}_2\\\vdots \\ \bar{T}^N\end{matrix}  \right\}\\ \displaystyle {\bf N}_v = [{\bf N}_1, {\bf N}_2,\cdots , {\bf N}_N ]^T ~,~ {\bf N}_p = {\bf N}_T = [{N}_1, { N}_2,\cdots , {N}_N ] \end{array} </math>
415
|}
416
| style="width: 5px;text-align: right;" | (25)
417
|}
418
419
with <math display="inline">{\bf N}_j = N_j {\bf I}_n</math> where <math display="inline">{\bf I}_n</math> is the <math display="inline">n\times n</math> unit matrix.
420
421
In Eq.([[#eq-25|25]]) vectors <math display="inline">\bar{v}</math>, <math display="inline">\bar {p}</math> and <math display="inline">\bar {T}</math>  contain the nodal velocities, the nodal pressures and the nodal temperatures for the whole mesh, respectively and the upperindex denotes the nodal value for each vector or scalar magnitude.
422
423
Substituting Eqs.([[#eq-24|24]]) into Eqs.([[#eq-17|17]]), ([[#eq-22|22]]) and ([[#eq-25|25]]) and choosing a Galerkin formulation with <math display="inline">w_i = q =  \hat w_i =N_i</math> leads to the following system of algebraic equations
424
425
<span id="eq-26"></span>
426
<span id="eq-26a"></span>
427
{| class="formulaSCP" style="width: 100%; text-align: left;" 
428
|-
429
| 
430
{| style="text-align: left; margin:auto;" 
431
|-
432
| style="text-align: center;" | <math>{\bf M}_0 {\dot{\bar{\bf v}}} + {\bf K}\bar{\bf v}+{\bf Q}\bar {\bf p}- {\bf f}_v={\bf 0} </math>
433
|}
434
| style="width: 5px;text-align: right;" | (26a)
435
|}
436
437
<span id="eq-26b"></span>
438
{| class="formulaSCP" style="width: 100%; text-align: left;" 
439
|-
440
| 
441
{| style="text-align: left; margin:auto;" 
442
|-
443
| style="text-align: center;" | <math>\textbf{M}_1 {\dot{\bar{\bf p}}}+\textbf{M}_2 {\ddot{\bar{\bf p}}}-{\bf Q}^T \bar{\bf v} + ({\bf L}+{\bf M}_b)\bar {\bf p}- {\bf f}_p={\bf 0} </math>
444
|}
445
| style="width: 5px;text-align: right;" | (26b)
446
|}
447
448
<span id="eq-26c"></span>
449
{| class="formulaSCP" style="width: 100%; text-align: left;" 
450
|-
451
| 
452
{| style="text-align: left; margin:auto;" 
453
|-
454
| style="text-align: center;" | <math>{\bf C} {\dot{\bar{\bf T}}}+\hat {\bf L} \bar{\bf T} -  {\bf f}_T={\bf 0}  </math>
455
|}
456
| style="width: 5px;text-align: right;" | (26c)
457
|}
458
459
where <math display="inline">{\dot{\bar{\bf a}}}</math> and <math display="inline">{\ddot{\bar{\bf a}}}</math> denote the first and second material time derivatives of the components of  a vector <math display="inline">{a}</math>. The different matrices and vectors in Eqs.([[#eq-26|26]]) are assembled from the element contributions given in Box 1.
460
461
<br/>
462
463
'''Remark 2.''' The boundary terms of vector <math display="inline">{\bf f}_p</math> can be incorporated in the matrices of Eq.([[#eq-26b|26b]]). This leads to a non symmetrical set of equations. These boundary terms are computed here iteratively within the incremental solution scheme.
464
465
'''Remark 3.''' The presence of matrix <math display="inline">{\bf M}_b</math> in Eq.([[#eq-26b|26b]]) 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 leads to considerable mass losses in viscous flows <span id='citeF-15'></span>[[#cite-15|[15]]].
466
467
'''Remark 4.''' For ''transient problems'' the stabilization parameter <math display="inline">\tau </math> of Eq.([[#eq-14|14]]) is computed 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
468
469
<span id="eq-27"></span>
470
{| class="formulaSCP" style="width: 100%; text-align: left;" 
471
|-
472
| 
473
{| style="text-align: left; margin:auto;" 
474
|-
475
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{(l^e)^2} + \frac{2\rho }{\Delta t}  \right)^{-1} </math>
476
|}
477
| style="width: 5px;text-align: right;" | (27)
478
|}
479
480
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(\Omega ^e)^{1/n_s}</math> where <math display="inline">\Omega ^e</math> is the element area (for 3-noded triangles) or volume (for 4-noded tetrahedra). For fluids with heterogeneous material the values of <math display="inline">\mu </math> and <math display="inline">\rho </math>  are computed at the element center.
481
482
For ''steady state'' problems the stabilization parameter is computed with Eq.([[#eq-27|27]]) substituting <math display="inline">\Delta t</math> by <math display="inline">\frac{l^e}{\vert {\bf v}^e\vert }</math>.  The characteristic boundary length <math display="inline">h_n</math> in the expression of <math display="inline">{\bf f}_p</math> (Box 1) has been taken equal to <math display="inline">l^e</math> in our computations.
483
484
{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
485
|-
486
| style="border-left: 2px solid;border-right: 2px solid;border-top: 2px solid;" | <math>\begin{array}{l}  \displaystyle \mathbf{M}_{0_{ij}}^e  =\int _{\Omega ^e}\rho {N}_{i}^e {N}_{j} {\bf I}_3 d\Omega ~,~\mathbf{K}^e_{ij} =\int _{\Omega ^e} \mathbf{B}_i^{eT} \mathbf{D}\mathbf{B}_j^e d\Omega ~,~\mathbf{Q}^e_{ij}  =\int _{\Omega ^e} \mathbf{B}_i^{eT} \mathbf{m} {N}_j^e d\Omega \\[.4cm] \displaystyle {M}_{1_{ij}}^e =\int _{\Omega ^e} \frac{1}{\kappa }{N}_i^e {N}_j^e d\Omega ~,~ {M}_{2_{ij}}^e =\int _{\Omega ^e} \frac{\tau }{ c^2 }{N}_i^e{N}_j^e d\Omega ~,~ {M}_{b_{ij}}^e = \int _{\Gamma _t} \frac{2\tau }{h_n} {N}_i^e{N}_j^e  d\Gamma \\[.4cm] 
487
\displaystyle {L}^e_{ij}= \int _{\Omega ^e} \tau ({\boldsymbol \nabla }^T {N}_i^e) {\boldsymbol \nabla }  N_j^e d\Omega ~,~ \mathbf{f}^e_{v_i}=  \int _{\Omega ^e}\mathbf{N}_{i}^e \mathbf{b} d\Omega + \int _{\Gamma _t}  \mathbf{N}_{i}^e {\bf t}^p d\Gamma \\[.4cm] 
488
\displaystyle {f}^e_{p_i}=\int _{\Gamma _t}\tau N_i^e \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} (2\mu \varepsilon _n -t_n)\right]d\Gamma - \int _{\Omega ^e} \tau  {\boldsymbol \nabla }^T {N}_i^e {b} d\Omega \\[.4cm] 
489
\displaystyle {C}^e_{ij}= \int _{\Omega ^e} \rho c N_i^e N_j^e d\Omega ~,~ \hat L^e_{ij}= \int _{\Omega ^e} k ({\boldsymbol \nabla }^T {N}_i^e) {\boldsymbol \nabla } N_j^e d\Omega ~,~ \displaystyle {f}^e_{T_i}= \int _{\Omega ^e} N_i^e Q d\Omega - \int _{\Gamma _q} N_i^e q^p_n d\Gamma \\[.4cm] \hbox{with } i,j=1,n.\\ \hbox{For 3D problems}\\[.4cm] 
490
\displaystyle{\bf B}_i^e = \left[\begin{matrix}  \displaystyle {\partial N_i^e \over \partial x_1} &0&0\\ \displaystyle{0}& \displaystyle {\partial N_i^e \over \partial x_2}&0\\ \displaystyle{0}&0&\displaystyle {\partial N_i^e \over \partial x_3}\\ \displaystyle {\partial N_i^e \over \partial x_2}&\displaystyle {\partial N_i^e \over \partial x_1}&0\\[.25cm] \displaystyle {\partial N_i^e \over \partial x_3}&0&\displaystyle {\partial N_i^e \over \partial x_1}\\[.25cm] \displaystyle{0}&\displaystyle {\partial N_i^e \over \partial x_3}&\displaystyle {\partial N_i^e \over \partial x_2}          \end{matrix}  \right]~,~ \mathbf{\bf N}_{i}^e = N_i^e {\bf I}_3 \quad \hbox{and} \quad  {\boldsymbol \nabla } = \left\{\begin{matrix}\displaystyle {\partial  \over \partial x_1} \\ \displaystyle {\partial  \over \partial x_2}\\\displaystyle {\partial  \over \partial x_3}  \end{matrix}  \right\} \end{array}</math>
491
|- style="border-bottom: 2px solid;"
492
| style="border-left: 2px solid;border-right: 2px solid;" | <math>N_i^e</math>: Local shape function of node <math display="inline">i</math> of element <math display="inline">e</math> [34,44]
493
|}
494
495
<div class="center" style="font-size: 75%;">'''Box 1'''. Element form of the matrices and vectors in Eqs.(26)</div>
496
497
==6 TRANSIENT SOLUTION OF THE DISCRETIZED EQUATIONS==
498
499
Eqs.([[#eq-26|26]]) are solved in time with an implicit Newton-Raphson type iterative scheme <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-44'></span><span id='citeF-46'></span>[[#cite-3|[3,7,44,46]]]. The basic steps within a time increment <math display="inline">[n,n+1]</math> are:<br/>
500
501
* Initialize variables: <math display="inline">({}^{n+1}{\bf x}^1,{}^{n+1}\bar {\bf v}^1,{}^{n+1}\bar{\bf p}^1,  {}^{n+1}\bar{\bf T}^1, {}^{n+1}\bar{\bf r}^1_m)\equiv  \left\{{}^{n}{\bf x},{}^{n}\bar{\bf v},{}^{n}\bar{\bf p},{}^{n}\bar{\bf T}, {}^{n}\bar{\bf r}_m \right\}</math>.
502
* Iteration loop: <math display="inline">i=1,NITER</math>.<p> For each iteration. 
503
504
</p>
505
506
===''Step 1. Compute the nodal velocity increments <math>\Delta \bar {v}</math>''===
507
508
From Eq.([[#eq-26a|26a]]), we deduce
509
510
<span id="eq-28a"></span>
511
{| class="formulaSCP" style="width: 100%; text-align: left;" 
512
|-
513
| 
514
{| style="text-align: left; margin:auto;" 
515
|-
516
| style="text-align: center;" | <math>{}^{n+1}{\bf H}_v^i \Delta \bar {\bf v} = - {}^{n+1}\bar{{\bf r}}_m^i \rightarrow  \Delta \bar{\bf v} </math>
517
|}
518
| style="width: 5px;text-align: right;" | (28a)
519
|}
520
521
with the momentum residual <math display="inline">\bar{\bf r}_m</math> and the iteration matrix <math display="inline">{\bf H}_v</math>  given by
522
523
<span id="eq-28b"></span>
524
{| class="formulaSCP" style="width: 100%; text-align: left;" 
525
|-
526
| 
527
{| style="text-align: left; margin:auto;" 
528
|-
529
| style="text-align: center;" | <math>\bar{\textbf{r}}_m = \textbf{M}_0 {\dot{\bar{\bf v}}} + \textbf{K}\bar{\bf v}+\textbf{Q}\bar {\bf p}-\textbf{f}_v\quad , \quad \textbf{H}_v = \frac{1}{\Delta t} \textbf{M}_0 + \textbf{K} + \textbf{K}_v  </math>
530
|}
531
| style="width: 5px;text-align: right;" | (28b)
532
|}
533
534
===''Step 2. Update the nodal velocities''===
535
536
<span id="eq-29"></span>
537
{| class="formulaSCP" style="width: 100%; text-align: left;" 
538
|-
539
| 
540
{| style="text-align: left; margin:auto;" 
541
|-
542
| style="text-align: center;" | <math>{}^{n+1}\bar{\bf v}^{i+1}= {}^{n+1}\bar{\bf v}^i + \Delta \bar{\bf v} </math>
543
|}
544
| style="width: 5px;text-align: right;" | (29)
545
|}
546
547
===''Step 3. Compute the nodal pressures'' <math> {}^{n+1}\bar{\bf p}^{i+1}</math>===
548
549
From Eq.([[#eq-26b|26b]]) we obtain
550
551
<span id="eq-30a"></span>
552
{| class="formulaSCP" style="width: 100%; text-align: left;" 
553
|-
554
| 
555
{| style="text-align: left; margin:auto;" 
556
|-
557
| style="text-align: center;" | <math>{}^{n+1} \textbf{H}_p^i  {}^{n+1}\bar{\bf p}^{i+1}= \frac{1}{\Delta t} \textbf{M}_1 {}^{n+1}\bar{\bf p}^i +  \frac{1}{\Delta t^2} \textbf{M}_2 (2 ^n \bar{\bf p} - ^{n-1} \bar{\bf p}) + {\bf Q}^T {}^{n+1}\bar{\bf v}^{i+1}+ {}^{n+1}\bar{\bf f}^{i}_p  \rightarrow  {}^{n+1}\bar{\bf p}^{i+1}   </math>
558
|}
559
| style="width: 5px;text-align: right;" | (30a)
560
|}
561
562
with
563
564
<span id="eq-30b"></span>
565
{| class="formulaSCP" style="width: 100%; text-align: left;" 
566
|-
567
| 
568
{| style="text-align: left; margin:auto;" 
569
|-
570
| style="text-align: center;" | <math>\textbf{H}_p = \frac{1}{\Delta t} \textbf{M}_1+ \frac{1}{\Delta t^2} \textbf{M}_2+\textbf{L} + \textbf{M}_b </math>
571
|}
572
| style="width: 5px;text-align: right;" | (30b)
573
|}
574
575
===''Step 4. Update the nodal coordinates''===
576
577
<span id="eq-31"></span>
578
{| class="formulaSCP" style="width: 100%; text-align: left;" 
579
|-
580
| 
581
{| style="text-align: left; margin:auto;" 
582
|-
583
| style="text-align: center;" | <math>{}^{n+1}{\bf x}^{i+1}= {}^{n+1}{\bf x}^i + \frac{1}{2} ({}^{n+1}\bar{\bf v}^{i+1} + {}^{n}\bar{\bf v})\Delta t </math>
584
|}
585
| style="width: 5px;text-align: right;" | (31)
586
|}
587
588
A more accurate expression for computing <math display="inline">{}^{n+1}{\bf x}^{i+1}</math> can be used involving the  nodal accelerations <span id='citeF-40'></span>[[#cite-40|[40]]].
589
590
===''Step 5.Compute the nodal temperatures''===
591
592
From Eq.([[#eq-26c|26c]]) we obtain
593
594
<span id="eq-32"></span>
595
{| class="formulaSCP" style="width: 100%; text-align: left;" 
596
|-
597
| 
598
{| style="text-align: left; margin:auto;" 
599
|-
600
| style="text-align: center;" | <math>\left[\frac{1}{\Delta t}{\bf C} + \hat{\bf L}   \right]\Delta \bar{\bf T}=- {}^{n+1} \bar{\bf r}^i_T \quad ,\quad  {}^{n+1}\bar{\bf T}^{i+1}= {}^{n}\bar{\bf T}^i + \Delta \bar{\bf T} </math>
601
|}
602
| style="width: 5px;text-align: right;" | (32)
603
|}
604
605
with
606
607
<span id="eq-33"></span>
608
{| class="formulaSCP" style="width: 100%; text-align: left;" 
609
|-
610
| 
611
{| style="text-align: left; margin:auto;" 
612
|-
613
| style="text-align: center;" | <math>{\bf r}_T= {\bf C} {\dot{\bar{\bf T}}} + \hat{\bf L} \bar{\bf T} - {\bf f}_T </math>
614
|}
615
| style="width: 5px;text-align: right;" | (33)
616
|}
617
618
===''Step 6. Check convergence''===
619
620
Verify the following conditions:
621
622
<span id="eq-34"></span>
623
{| class="formulaSCP" style="width: 100%; text-align: left;" 
624
|-
625
| 
626
{| style="text-align: left; margin:auto;" 
627
|-
628
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \Vert {}^{n+1}\bar {\bf v}^{i+1}- {}^{n+1}\bar{\bf v}^{i}\Vert \le e_v \Vert {}^{n}\bar {\bf v}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {\bf p}^{i+1}-{}^{n+1}\bar {\bf p}^i\Vert \le e_p  \Vert {}^{n}\bar{\bf p}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {\bf T}^{i+1} - {}^{n+1}\bar {\bf T}^i\Vert \le e_T  \Vert {}^{n}\bar{\bf T}\Vert  \end{array}   </math>
629
|}
630
| style="width: 5px;text-align: right;" | (34)
631
|}
632
633
where <math display="inline">e_v</math>, <math display="inline">e_p</math> and <math display="inline">e_T</math> are prescribed error norms. In our examples we have set <math display="inline">e_v = e_p=e_T = 10^{-3}</math>.
634
635
If  conditions ([[#eq-34|34]]) are satisfied then make <math display="inline">n \leftarrow n+1 </math> and proceed to the next  time step. Otherwise, make the iteration counter <math display="inline">i \leftarrow i+1 </math> and repeat Steps 1&#8211;5.
636
637
'''Remark 5.''' In Eqs.([[#eq-28|28]])&#8211;([[#eq-34|34]]) <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 this work the derivatives and integrals in all the matrices and the residual vectors <math display="inline">\bar{r}_m</math>  and <math display="inline">\bar{\bf r}_T</math> are computed on the ''discretized geometry at time'' <math display="inline">n</math>  while the nodal force vectors <math display="inline">{\bf f}_v</math>, <math display="inline">{\bf f}_p</math> and <math display="inline">{\bf f}_v</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-45'></span><span id='citeF-46'></span>[[#cite-3|[3,45,46]]].
638
639
'''Remark 6.'''  Including the  bulk stiffness matrix <math display="inline">\textbf{K}_v</math> in <math display="inline">\textbf{H}_v</math> has proven to be essential for the fast convergence, mass preservation and overall accuracy of the iterative solution <span id='citeF-10'></span><span id='citeF-41'></span>[[#cite-10|[10,41]]].     The element expression of <math display="inline">\textbf{K}_v</math> can be obtained as <span id='citeF-41'></span>[[#cite-41|[41]]]
640
641
<span id="eq-35"></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>{bf K}_v^e = \int _{\Omega ^e} \textbf{B}^T \textbf{m} \theta \Delta t\kappa \textbf{m}^T \textbf{B}d\Omega  </math>
648
|}
649
| style="width: 5px;text-align: right;" | (35)
650
|}
651
652
where <math display="inline">\theta </math> is a positive number such that <math display="inline">0<   \theta \le 1</math> that has the role of preventing the ill-conditioning of the iteration matrix <math display="inline">{\bf H}_v</math> for highly incompressible fluids. 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 <span id='citeF-10'></span>[[#cite-10|[10]]]. For fully incompressible fluids (<math display="inline">c</math> and <math display="inline">\kappa =\infty </math>), a finite  value of <math display="inline">\kappa </math> is used in practice in <math display="inline">{\bf K}_v</math> as this helps to obtaining an accurate solution for velocities and pressure  with reduced mass loss in few iterations per time step <span id='citeF-10'></span>[[#cite-10|[10]]]. These considerations, however, do not affect the value of <math display="inline">\kappa </math> within matrix <math display="inline">{\bf M}_1</math> in Eq.([[#eq-26b|26b]]) that vanishes for a fully incompressible fluid. Clearly, the value of the terms of <math display="inline">{\bf K}_v^e</math> can also be limited by reducing the time step size. This, however, leads to an increase in the cost of the computations. A similar approach for improving mass conservation in incompressible flows was proposed  in <span id='citeF-42'></span>[[#cite-42|[42]]].
653
654
'''Remark 7.''' The iteration matrix <math display="inline">\textbf{H}_v</math> in Eq.([[#eq-28a|28a]]) is an ''approximation of the exact tangent matrix'' in the updated Lagrangian formulation for a quasi-incompressible fluid <span id='citeF-40'></span>[[#cite-40|[40]]]. The simplified form of <math display="inline">\textbf{H}_v</math>  used in this work has yielded good results  with convergence achieved for the nodal velocities,  the pressure and the temperature in 3&#8211;4 iterations in all the problems analyzed.
655
656
'''Remark 8.''' The time step within a time interval <math display="inline">[n,n+1]</math> is chosen as <math display="inline">\Delta t =\min \left(\frac{^n l_{\min }^e}{\vert {}^n{\bf 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 Remark 4, <math display="inline">\vert ^n{\bf 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 defined as <math display="inline">\Delta t_b = \min \left(\frac{^n l_b}{\vert ^n{\bf v}_b\vert _{\max }}\right)</math> where <math display="inline">^n l_b</math> is the distance from the node to the boundary and <math display="inline">^n{\bf v}_b</math> is the velocity of the node. This definition of <math display="inline">\Delta t</math> intends that no node crosses a solid boundary during a time step.
657
658
A method that allows using large time steps in the integration of the PFEM equations can be found in <span id='citeF-16'></span>[[#cite-16|[16]]].
659
660
==7 ABOUT THE PARTICLE FINITE ELEMENT METHOD (PFEM)==
661
662
===7.1 The basis of the PFEM===
663
664
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 ''particles''. The 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-45'></span>[[#cite-3|[3,45]]].
665
666
The solution steps within a time step in the PFEM are as follows:
667
668
<div id='img-1'></div>
669
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
670
|-
671
|[[Image:draft_Samper_264539605-Figure2a.png|600px|Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time n   (t=ⁿt)  to   time n+2 (t=ⁿt +2∆t)  ]]
672
|- style="text-align: center; font-size: 75%;"
673
| colspan="1" | '''Figure 1:''' Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time <math>n</math>   (<math>t=^n t</math>)  to   time <math>n+2</math> (<math>t=^n t +2\Delta t</math>)  
674
|}
675
676
<ol>
677
678
<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-1|1]]). </li>
679
680
<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-8'></span>[[#cite-8|[8]]] 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-12'></span><span id='citeF-28'></span>[[#cite-12|[12,28]]]. </li>
681
682
<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-11'></span><span id='citeF-12'></span>[[#cite-11|[11,12]]]. </li>
683
684
<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>
685
686
<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>
687
688
<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>
689
690
</ol>
691
692
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.
693
694
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, while the solution of the equations (with typically 3 iterations per time step) and the system assembly consume approximately 70% and 15% of the CPU time per time step, respectively. These figures refer to analyses in a single processor Pentium IV PC <span id='citeF-36'></span>[[#cite-36|[36]]]. Considerable speed can be gained using parallel computing techniques.
695
696
In this work we will apply the PFEM to problems involving a rigid domain containing fluid particles only. Application of the PFEM in fluid and solid mechanics and in fluid-structure interaction problems can be found in <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-42'></span><span id='citeF-43'></span>[[#cite-4|[4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,28,29,35,36,39,40,41,42,43]]], as well in www.cimne.com/pfem.
697
698
==8 EXAMPLES==
699
700
===8.1 Sloshing of water in prismatic tank===
701
702
The problem has been solved first in 2D. Figure [[#img-2|2]] shows the analysis data. The fluid oscillates due to the hydrostatic forces induced by its original position.
703
704
<div id='img-2'></div>
705
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
706
|-
707
|[[Image:draft_Samper_264539605-Figure2_SL2D_official_input.png|570px|2D analysis of sloshing of water in rectangular tank. Initial geometry, analysis data and mesh of 5064 3-noded triangles discretizing the water in the tank]]
708
|- style="text-align: center; font-size: 75%;"
709
| colspan="1" | '''Figure 2:''' 2D analysis of sloshing of water in rectangular tank. Initial geometry, analysis data and mesh of 5064 3-noded triangles discretizing the water in the tank
710
|}
711
712
The problem has been run using different values of the parameter <math display="inline">\theta </math> in the tangent bulk stiffness matrix <math display="inline">{\bf K}_v^e</math> (Eq.([[#eq-41|41]])). The first set of results (Figures [[#img-3|3]] and [[#img-4|4]]) were obtained with <math display="inline">\theta =1</math>. The problem was then solved for <math display="inline">\theta = 0.08</math>, thereby, reducing in one order the magnitude the diagonal terms in <math display="inline">{\bf K}_v^e</math>.
713
714
Figure  [[#img-3|3]] shows snapshots of the water geometry at different times. Pressure contours are superposed to the deformed geometry of the fluid in the figures.
715
716
Figure [[#img-4|4]] shows the evolution of the percentage of water volume (i.e. mass) loss introduced by the numerical solution scheme. The accumulated volume loss (in percentage versus the initial volume) for the method proposed with <math display="inline">\theta =1</math> is approximately  1.33% over 20 seconds of simulation time (Figure [[#img-4|4]]a). The average volume variation in absolute value per time step is <math display="inline">1.09 \times 10^{-4}%</math> (Figure [[#img-4|4]]b). The total water volume loss is the sum of the losses induced by the numerical scheme and the losses due to the updating of the free surface  using the PFEM. No correction of mass was introduced at the end of each time step. Taking all this into account, the fluid volume loss over the analysis period is remarkably low.
717
718
The volume losses induced by the free surface updating can be reduced using a finer mesh in that region in conjunction with an enhanced alpha shape technique.
719
720
The total fluid volume loss can be reduced to almost zero by introducing a small correction in the free surface at the end of each time step <span id='citeF-41'></span>[[#cite-41|[41]]].
721
722
The fluid volume losses obtained using a standard first order fractional step method <span id='citeF-41'></span>[[#cite-41|[41]]]  and the PFEM are shown in Figure [[#img-4|4]]a  for comparison. Clearly the method proposed in this work leads to a reduction in the overall fluid volume loss, as well as in the volume loss per time step.
723
724
<div id='img-3'></div>
725
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
726
|-
727
|[[Image:draft_Samper_264539605-Figure4a_SL2D_t=5-7s.png|400px|]]
728
|[[Image:draft_Samper_264539605-Figure4b_SL2D_t=7-4s.png|400px|]]
729
|-
730
| style="text-align: center;font-size: 75%;"|(a) t=5.7s
731
| style="text-align: center;font-size: 75%;"| (b) t=7.4s
732
|-
733
|[[Image:draft_Samper_264539605-Figure4c_SL2D_t=13-3s.png|400px|]]
734
|[[Image:draft_Samper_264539605-Figure4d_SL2D_t=18-6s.png|400px|]]
735
|-
736
| style="text-align: center;font-size: 75%;"|(c) t=13.3s
737
| style="text-align: center;font-size: 75%;padding-bottom:10px;"| (d) t=18.6s
738
|- style="text-align: center; font-size: 75%;"
739
| colspan="2" | '''Figure 3:''' 2D sloshing of water in rectangular tank. Snapshots of water geometry at two different times (<math>\theta = 1</math>). Colours indicate pressure contours
740
|}
741
742
<div id='img-4'></div>
743
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
744
|-
745
|[[Image:draft_Samper_264539605-Figure5a_SL2D_allLosses.png|600px|]]
746
|-
747
| style="text-align: center;font-size: 75%;"|(a) Accumulated volume loss over 20 seconds of analysis
748
|-
749
|[[Image:draft_Samper_264539605-Figure5b_SL2D_lossesPerStep.png|600px|]]
750
|-
751
| style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Volume variation (in %) per time step over 20 seconds of analysis (Current method with <math>\theta = 1</math>)
752
|- style="text-align: center; font-size: 75%;"
753
| colspan="1" | '''Figure 4:''' 2D sloshing of water in rectangular tank. (a) Time evolution of the percentage of water volume loss  due to the numerical algorithm. (b) Average volume variation per time step. Current method. 1.09<math>\times 10^{-4}</math>% Fractional step: 2.07<math>\times 10^{-4}</math>% 
754
|}
755
756
757
Figure [[#img-5|5]] shows a comparison between the fluid volume loss for <math display="inline">\theta =1</math> and <math display="inline">\theta = 0.08</math> using the same time step in both cases (<math display="inline">\Delta t = 10^{-3}s</math>). Results  show that the reduction of the tangent bulk stiffness matrix terms leads to an improvement in the preservation of the initial volume of the fluid. It is noted that the convergence of the iterative solution for <math display="inline">\theta = 0.08</math> was the same as for <math display="inline">\theta =1</math>.
758
759
Figure [[#img-6|6]] shows that a similar improvement in the volume preservation can be obtained using <math display="inline">\theta = 1</math> and reducing the time step to <math display="inline">\Delta t = 10^{-4}s</math>. This, however, increases the cost of the computations.
760
761
These results indicate that accurate numerical results with reduced volume losses can be obtained by appropriately adjusting the parameter <math display="inline">\theta </math> in the tangent bulk modulus matrix while keeping the time step size to competitive values in terms of CPU cost. A study of the influence of <math display="inline">\theta </math> in the numerical solution for quasi-incompressible free surface fluids in terms of volume preservation and overall accuracy using the formulation here presented can be found in <span id='citeF-10'></span>[[#cite-10|[10]]].
762
763
More results for this example can be found in <span id='citeF-41'></span>[[#cite-41|[41]]].
764
765
<div id='img-5'></div>
766
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
767
|-
768
|[[Image:draft_Samper_264539605-Comparison-value.png|600px|2D sloshing of water in rectangular tank. Time evolution of percentage of water volume loss obtained using the current method with θ=0.08 (curve A) and θ= 1 (curve B) ∆t = 10⁻³s]]
769
|- style="text-align: center; font-size: 75%;"
770
| colspan="1" | '''Figure 5:''' 2D sloshing of water in rectangular tank. Time evolution of percentage of water volume loss obtained using the current method with <math>\theta =0.08</math> (curve A) and <math>\theta = 1</math> (curve B) <math>\Delta t = 10^{-3}s</math>
771
|}
772
773
<div id='img-6'></div>
774
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
775
|-
776
|[[Image:draft_Samper_264539605-volume-preservation.png|600px|2D sloshing of water in rectangular tank. Time evolution of percentage of water volume loss obtained with the current method. Curve A: θ=1 and ∆t = 10⁻⁴s. Curve B: θ= 1 and ∆t = 10⁻³s]]
777
|- style="text-align: center; font-size: 75%;"
778
| colspan="1" | '''Figure 6:''' 2D sloshing of water in rectangular tank. Time evolution of percentage of water volume loss obtained with the current method. Curve A: <math>\theta =1</math> and <math>\Delta t = 10^{-4}s</math>. Curve B: <math>\theta = 1</math> and <math>\Delta t = 10^{-3}s</math>
779
|}
780
781
782
Figures [[#img-7|7]] and [[#img-8|8]] show a similar set of results for the 3D analysis of the same sloshing problem using a relative coarse initial mesh of 106771 4-noded tetrahedra and <math display="inline">\theta = 1</math>. It is remarkable that the percentage of total fluid volume loss due to the numerical scheme after 10 seconds of analysis is approximately 1%.
783
784
785
<div id='img-7'></div>
786
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
787
|-
788
|colspan="2" | [[File:Draft_Samper_264539605_2248_Figure7a_SL3D_AS2D_input.png]]
789
|150px|]]
790
|-
791
|[[Image:draft_Samper_264539605-Figure7bSL3D_AS2D_5_7.png|300px|]]
792
|[[Image:draft_Samper_264539605-Figure7c_SL3D_AS2D_7_4.png|300px|]]
793
|-
794
| style="text-align: center;font-size: 75%;"|(a) t=5.7s
795
| style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) t=7.4s
796
|- style="text-align: center; font-size: 75%;"
797
| colspan="2" | '''Figure 7:''' 3D analysis of sloshing of water in prismatic tank (<math>\theta = 1</math>). Analysis data and snapshots of water geometry at <math>t=5.7</math>s (a) and <math>t=7.4</math>s (b)
798
|}
799
800
<div id='img-8'></div>
801
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
802
|-
803
|[[Image:draft_Samper_264539605-Fig13b_New.png|500px]]
804
|-
805
|style="text-align: center;font-size: 75%;"| (a)
806
|-
807
|[[Image:draft_Samper_264539605-Fig13c_New.png|500px]]
808
|-
809
|style="text-align: center;font-size: 75%;padding-bottom:10px;"| (b) 
810
|- style="text-align: center; font-size: 75%;"
811
| colspan="1" | '''Figure 8:''' 3D analysis of sloshing of water in prismatic tank (<math>\theta = 1</math>). (a)  Time evolution of accumulated water volume loss (in percentage) due to the numerical algorithm. (b) Volume loss (in %) per time step over 2 seconds of analysis. Average volume loss per time step: 1.64<math>\times 10^{-4}</math>%
812
|}
813
814
815
===8.2 Falling of a water sphere in a cylindrical tank containing water===
816
817
This example is the 3D analysis of the impact of a sphere made of water as it falls in a cylindrical tank containing water. Both the water in the sphere and in the tank mix in a single fluid after the impact. Figure [[#img-9|9]] shows the material and analysis data and the initial discretization of the sphere and the water in the tank in 88892 4-noded tetrahedra. The problem was solved with the new stabilized method presented in the paper with <math display="inline">\theta = 1</math>. Figure [[#img-10|10]] shows snapshots of the mixing process at different times. An average of four  iterations for convergence of the velocity and the pressure were needed during all the steps of the analysis. The total water mass lost in the sphere and the tank due to the numerical algorithm was <math display="inline">\simeq </math> 2% after 3 seconds of analysis (Figure [[#img-11|11]]a).
818
819
820
<div id='img-9'></div>
821
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
822
|-
823
|[[Image:draft_Samper_264539605-Figure19_ball_3D_input.png|600px|Water sphere falling  in a tank filled with water. Analysis data and initial discretization of the sphere and the water in the tank with 88892 4-noded tetrahedra]]
824
|- style="text-align: center; font-size: 75%;"
825
| colspan="1" | '''Figure 9:''' Water sphere falling  in a tank filled with water. Analysis data and initial discretization of the sphere and the water in the tank with 88892 4-noded tetrahedra
826
|}
827
828
<div id='img-10'></div>
829
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
830
|-
831
|[[Image:draft_Samper_264539605-figure20_ball1.png|200px|]]
832
|[[Image:draft_Samper_264539605-figure20_ball2.png|200px|]]
833
|-
834
|style="text-align: center;font-size: 75%;"| (a) t=0.175s
835
|style="text-align: center;font-size: 75%;"| (b) t=0.275s
836
|-
837
|[[Image:draft_Samper_264539605-figure20_ball3.png|200px|]]
838
|[[Image:draft_Samper_264539605-figure20_ball4.png|200px|]]
839
|-
840
|style="text-align: center;font-size: 75%;"| (c) t=0.5s
841
|style="text-align: center;font-size: 75%;padding-buttom:10px;"| (d) t=0.9s
842
|- style="text-align: center; font-size: 75%;"
843
| colspan="2" | '''Figure 10:''' Water sphere falling  in tank containing water. Evolution of the impact and mixing of the two liquids at different times. Results for <math>\theta = 1</math>
844
|}
845
846
<div id='img-11'></div>
847
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
848
|-
849
|[[Image:draft_Samper_264539605-Fig26a_NEW.png|400px|]]
850
|-
851
|style="text-align: center;font-size: 75%;"| (a)
852
|-
853
|[[Image:draft_Samper_264539605-Fig26B_NEW.png|400px|]]
854
|-
855
|style="text-align: center;font-size: 75%;padding-bottom:10px;"| (b)
856
|- style="text-align: center; font-size: 75%;"
857
| colspan="1" | '''Figure 11:''' Water sphere falling in a tank containing water. (a) Accumulated volume over three seconds of analysis due to the numerical algorithm (<math>\theta =1</math>). (b) Volume loss (in %) per time step. Average volume variation per time step: 2.54<math>\times 10^{-4}</math>%
858
|}
859
860
861
===8.3 Sloshing of a fluid in a heated tank===
862
863
A fluid at initial temperature <math display="inline">T=20\,^{\circ } C</math> oscillates due to the hydrostatic forces induced by its initial position in a rectangular tank heated to a uniform and constant temperature of <math display="inline">T=75\,^{\circ } C</math>. The geometry and the problem data of the 2D simulation are shown in Figure [[#img-12|12]]. The fluid, with a very high thermal conductivity, changes its temperature only due to the contact with the hotter tank walls. The heat flux along the free surface has been considered null. The fluid domain has been initially discretized with 2828 3-noded triangles. The coupled thermal-fluid dynamics simulation has been run for 100 s using a time step increment of <math display="inline">\Delta t</math>=<math display="inline">0.005 s</math>.
864
865
Figure [[#img-13|13]] shows some snapshots of the numerical simulation. The temperature contours have been superposed on the fluid domain at the different time instants.
866
867
In Figure [[#img-14|14]] the evolution of temperature with time at the points <math display="inline">A, B</math> and <math display="inline">C</math> of  Figure [[#img-12|12]] is plotted. The coordinates of these sample points are (<math display="inline">m,0.1m</math>), (<math display="inline">m,0.4m</math>) and (<math display="inline">m,0.1m</math>), respectively. Figures [[#img-13|13]]  and [[#img-14|14]]  show that the fluid does not heat uniformly because of the convection effect automatically captured by the Lagrangian technique here presented.
868
869
<div id='img-12'></div>
870
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
871
|-
872
|[[Image:draft_Samper_264539605-ThermalSloshingInput.png|300px|]]
873
|[[Image:draft_Samper_264539605-ThermalSLData.png|200px|2D sloshing of a fluid in a heated tank. Initial geometry, problem data, thermal boundary and initial conditions]]
874
|- style="text-align: center; font-size: 75%;"
875
| colspan="2" | '''Figure 12:''' 2D sloshing of a fluid in a heated tank. Initial geometry, problem data, thermal boundary and initial conditions
876
|}
877
878
<div id='img-13'></div>
879
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
880
|-
881
|[[Image:draft_Samper_264539605-ThermalSloshing.png|500px|2D sloshing of a fluid in a heated tank. Snapshots of fluid geometry at six different times. Colours indicate temperature contours]]
882
|- style="text-align: center; font-size: 75%;"
883
| colspan="1" | '''Figure 13:''' 2D sloshing of a fluid in a heated tank. Snapshots of fluid geometry at six different times. Colours indicate temperature contours
884
|}
885
886
<div id='img-14'></div>
887
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
888
|-
889
|[[Image:draft_Samper_264539605-ThermalSloshingGraph.png|380px|2D sloshing of a fluid in a heated tank. Evolution of temperature with time at the points A, B and C of  Figure [[#img-12|12]]]]
890
|- style="text-align: center; font-size: 75%;"
891
| colspan="1" | '''Figure 14:''' 2D sloshing of a fluid in a heated tank. Evolution of temperature with time at the points <math>A, B</math> and <math>C</math> of  Figure [[#img-12|12]]
892
|}
893
894
895
===8.4 Falling of a cylindrical  object in a heated tank filled with fluid===
896
897
An elastic object falls in a tank containing a fluid at rest. The tank walls are maintained at temperature <math display="inline">T</math>=<math display="inline">75\,^{\circ } C</math> during the whole analysis, while the fluid and the solid object have an initial temperature of <math display="inline">T=20\,^{\circ } C</math>.  The geometry and the problem data of the 2D simulation as well the thermal initial and boundary conditions, are shown in Figure [[#img-15|15]]. Both the fluid and the solid have a high thermal conductivity. The heat flux along the fluid and solid surfaces in contact with the air has been considered null. The fluid and the solid domains have been discretized with 1986 and 108 3-noded triangular finite elements, respectively. The duration of the simulation is 10 s and the time step increment chosen is <math display="inline">\Delta t</math>=<math display="inline">0.0005 s</math>.
898
899
Figure [[#img-16|16]] collects some representative snapshots of the numerical simulation with the temperature results plotted over the fluid and the solid domains.
900
901
The graph of Figure [[#img-17|17]] is the evolution of temperature at the central point of the solid object. As expected, its temperature tends to <math display="inline">T</math>=<math display="inline">75\,^{\circ } C</math>.
902
903
<div id='img-15'></div>
904
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
905
|-
906
|[[Image:draft_Samper_264539605-ThermalBallInput.png|470px|]]
907
|-
908
|[[Image:draft_Samper_264539605-FSIData.png|500px|Falling of a solid object in a heated tank filled with fluid. Initial geometry, problem data, thermal boundary and initial conditions]]
909
|- style="text-align: center; font-size: 75%;"
910
| colspan="1" | '''Figure 15:''' Falling of a solid object in a heated tank filled with fluid. Initial geometry, problem data, thermal boundary and initial conditions
911
|}
912
913
<div id='img-16'></div>
914
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
915
|-
916
|[[Image:draft_Samper_264539605-ThermalBall.png|400px|Falling of a solid object in a heated tank filled with fluid. Snapshots at six different times. Colours indicate temperature contours]]
917
|- style="text-align: center; font-size: 75%;"
918
| colspan="1" | '''Figure 16:''' Falling of a solid object in a heated tank filled with fluid. Snapshots at six different times. Colours indicate temperature contours
919
|}
920
921
<div id='img-17'></div>
922
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
923
|-
924
|[[Image:draft_Samper_264539605-FSIGraph.png|400px|Falling of a solid object in a heated tank filled with fluid. Time evolution of the temperature  at the center of the solid.]]
925
|- style="text-align: center; font-size: 75%;"
926
| colspan="1" | '''Figure 17:''' Falling of a solid object in a heated tank filled with fluid. Time evolution of the temperature  at the center of the solid.
927
|}
928
929
==9 CONCLUDING REMARKS==
930
931
We have presented a new FIC-based stabilized Lagrangian finite element method for thermal-mechanical analysis of quasi and fully incompressible flows and FSI problems that has excellent mass preservation properties. The method has been successfully applied to the adiabatic and thermal-mechanical analysis of free-surface quasi-incompressible flows using the PFEM and an updated Lagrangian formulation. These problems are more demanding in terms of the mass preservation features of the numerical algorithm. The method proposed has yielded excellent results for 2D and 3D adiabatic and thermally-coupled free surface flow problems involving surface waves, water splashing, violent impact of flows with containment walls  and FSI situations.
932
933
==ACKNOWLEDGEMENTS==
934
935
This research was partially supported by the Advanced Grant  project SAFECON of the European Research Council.
936
937
==Bibliography==
938
939
<div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">
940
941
<div id="cite-1"></div>
942
[[#citeF-1|[1]]]  Aubry R, Idelsohn SR, Oñate E (2005)  Particle finite element method in fluid-mechanics including thermal convection-diffusion. Computers and Structures,  Vol. 83, (17-18), pp 1459&#8211;1475.
943
944
<div id="cite-2"></div>
945
[[#citeF-2|[2]]]   Aubry R, Idelsohn SR, Oñate E (2006)  Fractional step like schemes for free surface problems with thermal coupling using the Lagrangian PFEM. Computational Mechanics,  Vol. 38 (4-5), pp. 294&#8211;309.
946
947
<div id="cite-3"></div>
948
[[#citeF-3|[3]]]  Belytschko T, Liu WK, Moran B, Non linear finite element for continua and structures, 2d Edition, Wiley, 2013.
949
950
<div id="cite-4"></div>
951
[[#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
952
953
<div id="cite-5"></div>
954
[[#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
955
956
<div id="cite-6"></div>
957
[[#citeF-6|[6]]]  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
958
959
<div id="cite-7"></div>
960
[[#citeF-7|[7]]]   Donea J,  Huerta A  (2003) Finite element method for flow problems. J. Wiley
961
962
<div id="cite-8"></div>
963
[[#citeF-8|[8]]]  Edelsbrunner H,  Mucke EP (1999) Three dimensional alpha shapes. ACM   Trans. Graphics  13:43&#8211;72
964
965
<div id="cite-9"></div>
966
[[#citeF-9|[9]]]  Felippa F, Oñate E (2007)  Nodally exact Ritz discretizations of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Comput. Mech. 39:91&#8211;111
967
968
<div id="cite-10"></div>
969
[[#citeF-10|[10]]]  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. Submitted to Int. J. Numer. Meth. Biomed. Engng.
970
971
<div id="cite-11"></div>
972
[[#citeF-11|[11]]]  Idelsohn SR, Calvo N, Oñate E (2003c) Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng. 192(22-24):2649&#8211;2668
973
974
<div id="cite-12"></div>
975
[[#citeF-12|[12]]]  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. J. Numer. Meth. Biomed. Engng.,  61(7):964&#8211;989
976
977
<div id="cite-13"></div>
978
[[#citeF-13|[13]]]  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
979
980
<div id="cite-14"></div>
981
[[#citeF-14|[14]]]  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
982
983
<div id="cite-15"></div>
984
[[#citeF-15|[15]]]  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
985
986
<div id="cite-16"></div>
987
[[#citeF-16|[16]]]  Idelsohn SR, Nigro N, Limache A, Oñate E  (2012) Large time-step explicit integration method for solving problem with dominant convection. Comput. Methods Appl. Mech. Engrg. 217-220:168–-185
988
989
<div id="cite-17"></div>
990
[[#citeF-17|[17]]]  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
991
992
<div id="cite-18"></div>
993
[[#citeF-18|[18]]]  Limache A, Idelsohn, SR, Rossi R, Oñate E (2007) The violation of objectivity in Laplace formulation of the Navier-Stokes equations. Int. J. Numerical Methods in Fluids, 54:639&#8211;664.
994
995
<div id="cite-19"></div>
996
[[#citeF-19|[19]]]  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
997
998
<div id="cite-20"></div>
999
[[#citeF-20|[20]]]  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
1000
1001
<div id="cite-21"></div>
1002
[[#citeF-21|[21]]]  Oñate E, Manzan M (1999)  A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. International Journal for Numerical Methods in Fluids,  Vol. 31, pp. 203&#8211;221.
1003
1004
<div id="cite-22"></div>
1005
[[#citeF-22|[22]]]  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
1006
1007
<div id="cite-23"></div>
1008
[[#citeF-23|[23]]]   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
1009
1010
<div id="cite-24"></div>
1011
[[#citeF-24|[24]]]  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
1012
1013
<div id="cite-25"></div>
1014
[[#citeF-25|[25]]]  Oñate E, Taylor RL,  Zienkiewicz OC, Rojek J (2003) A residual correction method based on finite calculus. Engineering Computations 20:629&#8211;658
1015
1016
<div id="cite-26"></div>
1017
[[#citeF-26|[26]]]   Oñate E (2004) Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng. 60(1):255&#8211;281
1018
1019
<div id="cite-27"></div>
1020
[[#citeF-27|[27]]]   Oñate E, Rojek J, Taylor R, Zienkiewicz O (2004a) Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. Int. J. Numer. Meth. Engng. 59(11):1473&#8211;1500
1021
1022
<div id="cite-28"></div>
1023
[[#citeF-28|[28]]]  Oñate E, Idelsohn SR, Del Pin F, Aubry R (2004b) The particle   finite element method. An overview. Int. J. Comput. Methods    1(2):267&#8211;307
1024
1025
<div id="cite-29"></div>
1026
[[#citeF-29|[29]]]  Oñate E, M.A. Celigueta, Idelsohn SR (2006a) Modeling bed erosion in   free surface flows by the Particle Finite Element Method, Acta   Geotechnia 1(4):237&#8211;252
1027
1028
<div id="cite-30"></div>
1029
[[#citeF-30|[30]]]  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
1030
1031
<div id="cite-31"></div>
1032
[[#citeF-31|[31]]]  Oñate E, García J,  SR Idelsohn, F. Del Pin (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
1033
1034
<div id="cite-32"></div>
1035
[[#citeF-32|[32]]]  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
1036
1037
<div id="cite-33"></div>
1038
[33]  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
1039
1040
<div id="cite-34"></div>
1041
[[#citeF-34|[34]]]  Oñate E (2009)  Structural analysis with the finite element method. Linear statics. Volume 1. Basis and Solids.  CIMNE-Springer
1042
1043
<div id="cite-35"></div>
1044
[[#citeF-35|[35]]]  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
1045
1046
<div id="cite-36"></div>
1047
[[#citeF-36|[36]]]  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.
1048
1049
<div id="cite-37"></div>
1050
[[#citeF-37|[37]]]  Oñate E, Nadukandi P, Idelsohn SR, García J, Felippa C (2011) A family of residual-based stabilized finite element methods for Stokes flows. Int. J. Num. Methods in Fluids, 65 (1-3): 106&#8211;134
1051
1052
<div id="cite-38"></div>
1053
[[#citeF-38|[38]]]  Oñate E, Idelsohn SR, Felippa C (2011) Consistent pressure Laplacian stabilization for incompressible continua via higher-order finite calculus. Int. J. Numer. Meth. Engng, 87 (1-5): 171&#8211;195
1054
1055
<div id="cite-39"></div>
1056
[[#citeF-39|[39]]]  Oñate E, Nadukandi P, Idelsohn SR (2014) P1/P0+ elements for incompressible  flows  with discontinuous material properties. Comput. Meth. Appl. Mech. Engng. 271:185-209
1057
1058
<div id="cite-40"></div>
1059
[[#citeF-40|[40]]]  Oñate E, Carbonell JM  (2013) Updated Lagrangian finite element formulation for quasi-incompressible fluids. Research Report PI393, CIMNE. Submitted to Comput. Mechanics
1060
1061
<div id="cite-41"></div>
1062
[[#citeF-41|[41]]]  Oñate E, Franci A, Carbonell JM (2013) Lagrangian formulation  for finite element analysis of quasi-incompressible fluids with reduced mass losses. Int. J. Num. Meth. Fluids, DOI: 10.1002/fld.3870
1063
1064
<div id="cite-42"></div>
1065
[[#citeF-42|[42]]]  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
1066
1067
<div id="cite-43"></div>
1068
[[#citeF-43|[43]]]  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
1069
1070
<div id="cite-44"></div>
1071
[[#citeF-44|[44]]]  Zienkiewicz OC,  Taylor RL, Zhu JZ (2005) The finite element method. The basis,  6th Ed., Elsevier
1072
1073
<div id="cite-45"></div>
1074
[[#citeF-45|[45]]]  Zienkiewicz OC,  Taylor RL (2005) The finite element method for   solid and structural mechanics,  6th Ed., Elsevier
1075
1076
<div id="cite-46"></div>
1077
[[#citeF-46|[46]]]  Zienkiewicz OC,  Taylor RL, Nithiarasu P (2005)  The finite element   method for fluid dynamics,  6th Ed.,  Elsevier
1078
</div>
1079

Return to Onate et al 2014e.

Back to Top

Document information

Published on 01/01/2014

DOI: 10.1007/978-3-319-06136-8
Licence: CC BY-NC-SA license

Document Score

0

Views 83
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?