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 ''Comput. Meth. Appl. Mech. Engng.'', Vol. 293, pp. 191-206, 2015
2
3
4
==Abstract==
5
6
We present a 3-noded triangle and a 4-noded tetrahedra with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown ''constant pressure field'' and ''a prescribed linear field'' that satisfies the steady state momentum equations for a constant body force. The elements are termed P1/P0+ as the “effective” pressure field is linear, although the unknown pressure field is piecewise constant within each element. The elements  have an excellent behaviour for incompressible viscous flow problems with discontinuous material properties formulated in either Eulerian or Lagrangian descriptions. The necessary numerical stabilization for dealing with the inf-sup condition imposed by the incompressibility constraint and high convective effects (in Eulerian flows) is introduced via the Finite Calculus (FIC) approach. For the sake of clarity, the element derivation is presented first for the simpler Stokes equations written in the standard Eulerian frame. The extension of the formulation to the Navier-Stokes equations written in the Eulerian and Lagrangian frameworks is straightforward and is presented in the second part of the paper.
7
8
The efficiency and accuracy of the new P1/P0+ triangle is verified by solving a set of incompressible multifluid flow problems using a Lagrangian approach and a classical Eulerian description. The excellent performance of the new triangular element in terms of mass conservation and general accuracy for analysis of fluids with discontinuous material properties is highlighted.
9
10
'''Keywords''': P1/P0+ elements, Incompressible  flow,  Discontinuous material properties, Multifluids
11
12
13
==1 INTRODUCTION==
14
15
Preservation of mass is a great challenge in the numerical study of incompressible flow problems. Mass losses can be introduced by the so-called stabilization terms which are typically added to the discretized forms of the momentum and mass balance equations in order to account for large values of the convective acceleration terms in the momentum equations in high Reynolds number flows and to satisfy the inf-sup condition imposed by the incompressibility constraint when equal order interpolation of the velocities and the pressure are used in mixed finite element methods <span id='citeF-4'></span><span id='citeF-11'></span><span id='citeF-39'></span>[[#cite-4|[4,11,39]]].
16
17
Mass loss  also typically occurs in the numerical study of free-surface incompressible flow problems using Eulerian and Lagrangian descriptions. In both cases, the inaccuracy in capturing the free surface can lead to considerable mass losses unless special numerical schemes are used <span id='citeF-22'></span>[[#cite-22|[22]]].
18
19
An additional source of mass loss  occurs in the numerical analysis of the so-called multifluid problems with discontinuous changes in the viscosity and/or the density in parts of the domain. Most numerical methods have difficulties for accurately capturing  the jumps in the pressure and/or the pressure gradient at the interfaces between the different fluids <span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-25'></span>[[#cite-20|[20,21,25]]].
20
21
The numerical study of multifluids  via the FEM and similar computational techniques has been the subject of much research in last decades. Several authors have proposed alternative stabilized FEM procedures for accounting for the discontinuity in the pressure (and/or the pressure gradient) at the interface of fluids with different viscosity (and/or pressure). Among these formulations we note those based in injecting a discontinuous pressure field within the appropriate elements <span id='citeF-2'></span>[[#cite-2|[2]]] and those based on using a stabilized formulation based on the introduction of stabilization terms including the jumps in the pressure and the viscous terms at the element boundaries <span id='citeF-3'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-9'></span><span id='citeF-12'></span><span id='citeF-14'></span><span id='citeF-17'></span><span id='citeF-23'></span>[[#cite-3|[3,5,6,7,9,12,14,17,23]]].
22
23
In this work we present a new 3-noded triangle (and the 4-noded tetrahedron counterpart) with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown ''constant pressure field'' and ''a prescribed linear field'' that satisfies the steady state momentum equations for a constant body force. The so-called P1/P0+  elements  have an excellent behaviour for incompressible viscous flow problems formulated in Eulerian and Lagrangian descriptions. For the sake of clarity, the elements derivation is presented first for the simpler Stokes equations written in the standard Eulerian frame. The extension of the formulation to the Navier-Stokes equations written in the Eulerian and Lagrangian frameworks  is straightforward and is presented in the second part of the paper. A motivation of this work is the study of multifluids  problems using a Lagrangian formulation via the Particle Finite Element Method (PFEM) <span id='citeF-10'></span>[[#cite-10|[10]]],<span id='citeF-18'></span>[[#cite-18|[18]]]&#8211;<span id='citeF-21'></span>[[#cite-21|[21]]],<span id='citeF-25'></span><span id='citeF-32'></span>[[#cite-25|[25,32]]],<span id='citeF-36'></span>[[#cite-36|[36]]]&#8211;<span id='citeF-38'></span>[[#cite-38|[38]]] and similar procedures.
24
25
The success of the P1/P0+ formulation lays in the consistent derivation of a residual-based expression of the mass balance equation using the Finite Calculus (FIC) technique. The FIC approach in mechanics is based on expressing the equations of balance of mass and momentum in a domain of finite size and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain.  In addition to the standard terms of infinitesimal theory, the FIC forms of the balance equations contain derivatives of the classical differential equations in mechanics multiplied by characteristic distances in space and time. Examples of stabilized FIC-FEM formulations in fluid and solid mechanics can be found in <span id='citeF-26'></span>[[#cite-26|[26]]]&#8211;<span id='citeF-31'></span>[[#cite-31|[31]]],<span id='citeF-33'></span>[[#cite-33|[33]]]&#8211;<span id='citeF-35'></span>[[#cite-35|[35]]]. In our work we use the FIC forms of the mass balance equation in space and time for obtaining a variational residual form useful for finite element analysis.
26
27
The discretized stabilized variational form for the mass balance equation using a P1/P0+ approximation for the velocity/pressure variables involves the jumps in the viscous stresses, the pressure, the surface tractions and the acceleration term in the normal direction to each side (or face) of the elements. These stabilization terms resemble those proposed by other authors for similar fluid flow problems <span id='citeF-3'></span><span id='citeF-9'></span><span id='citeF-12'></span><span id='citeF-14'></span><span id='citeF-17'></span><span id='citeF-23'></span>[[#cite-3|[3,9,12,14,17,23]]]. The method presented in this paper yields a consistent and extended form of the stabilization terms that has proven to give a superior behaviour in terms of mass conservation and overall stability in multifluids problems.
28
29
The lay-out of the paper is the following. In the next section we present the basic equations for an incompressible Stokes fluid. Next we derive the stabilized variational FIC form of the mass balance equation. Then the P1/P0+ finite element discretization for the 3-noded triangle and the 4-noded tetrahedron is presented and the key matrices and vectors of the discretized system of  equations are given. Details of the solution of the FEM equations  are given.
30
31
The extension of the general stabilized FIC-FEM formulation to Navier-Stokes flows using a standard Eulerian approach and a Lagrangian description is outlined.
32
33
The efficiency and accuracy of the new P1/P0+ triangle is verified by solving a set of transient incompressible multifluid flow problems using a Lagrangian approach and a steady state problem via a classical Eulerian description. The excellent performance of the P1/P0+ triangle in terms of mass conservation and general accuracy for fluid flow problems with discontinuous material properties is highlighted.
34
35
==2 BASIC EQUATIONS==
36
37
We write the governing equations for an incompressible Stokes flow problem  as follows.
38
39
====''Momentum equations''====
40
41
<span id="eq-1"></span>
42
{| class="formulaSCP" style="width: 100%; text-align: left;" 
43
|-
44
| 
45
{| style="text-align: left; margin:auto;width: 100%;" 
46
|-
47
| style="text-align: center;" | <math>r_{m_i}:= \rho \frac{\partial v_i}{\partial t}-{\partial \sigma _{ij} \over \partial x_j}-b_i=0\quad , \quad i,j=1,n_s \quad \hbox{in } \Omega  </math>
48
|}
49
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
50
|}
51
52
In Eq.([[#eq-1|1]]), <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">n_s</math> is the number of space dimension (i.e. <math display="inline">n_s=3</math> for a 3D problem), <math display="inline">\Omega </math> is the analysis domain 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
53
54
<span id="eq-2"></span>
55
{| class="formulaSCP" style="width: 100%; text-align: left;" 
56
|-
57
| 
58
{| style="text-align: left; margin:auto;width: 100%;" 
59
|-
60
| style="text-align: center;" | <math>\sigma _{ij} = s_{ij}+ p \delta _{ij} </math>
61
|}
62
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
63
|}
64
65
where <math display="inline">\delta _{ij}</math> is the Kronecker delta. Note that the pressure is assumed to be positive for a tension state.
66
67
Summation of terms with repeated indices is assumed in Eq.([[#eq-1|1]]) and in  the following, unless otherwise specified.
68
69
====''Constitutive equation and volumetric strain rates''====
70
71
The relationship between the deviatoric stresses and the strain rates has the standard form for a Newtonian fluid,
72
73
<span id="eq-3"></span>
74
{| class="formulaSCP" style="width: 100%; text-align: left;" 
75
|-
76
| 
77
{| style="text-align: left; margin:auto;width: 100%;" 
78
|-
79
| style="text-align: center;" | <math>s_{ij} = 2\mu \left(\varepsilon _{ij} - {1 \over 3}\varepsilon _v \delta _{ij}\right) </math>
80
|}
81
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
82
|}
83
84
where  the strain rates <math display="inline">\varepsilon _{ij}</math> are related to the velocities by
85
86
<span id="eq-4"></span>
87
{| class="formulaSCP" style="width: 100%; text-align: left;" 
88
|-
89
| 
90
{| style="text-align: left; margin:auto;width: 100%;" 
91
|-
92
| 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>
93
|}
94
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
95
|}
96
97
In Eq.([[#eq-4|4]]) <math display="inline">\varepsilon _v</math> is the volumetric strain rate defined as
98
99
<span id="eq-5"></span>
100
{| class="formulaSCP" style="width: 100%; text-align: left;" 
101
|-
102
| 
103
{| style="text-align: left; margin:auto;width: 100%;" 
104
|-
105
| style="text-align: center;" | <math>\varepsilon _v= \varepsilon _{ii} </math>
106
|}
107
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
108
|}
109
110
Substituting Eqs.([[#eq-2|2]]) and ([[#eq-4|4]]) into ([[#eq-1|1]]) gives a  useful form of the momentum equations as
111
112
<span id="eq-6"></span>
113
{| class="formulaSCP" style="width: 100%; text-align: left;" 
114
|-
115
| 
116
{| style="text-align: left; margin:auto;width: 100%;" 
117
|-
118
| style="text-align: center;" | <math>\rho \frac{\partial v_i}{\partial t} - {\partial  \over \partial x_j}(2\mu \varepsilon _{ij})+{\partial  \over \partial x_i}\left(\frac{2}{3}\mu \varepsilon _v\right)- {\partial p \over \partial x_i}-b_i =0\quad , \quad i,j=1,n_s  </math>
119
|}
120
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
121
|}
122
123
====''Boundary conditions''====
124
125
The boundary conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann/traction (<math display="inline">\Gamma _t</math>) boundaries are
126
127
<span id="eq-7.a"></span>
128
{| class="formulaSCP" style="width: 100%; text-align: left;" 
129
|-
130
| 
131
{| style="text-align: left; margin:auto;width: 100%;" 
132
|-
133
| style="text-align: center;" | <math>v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v </math>
134
|}
135
| style="width: 5px;text-align: right;white-space: nowrap;" | (7.a)
136
|}
137
138
<span id="eq-7.b"></span>
139
{| class="formulaSCP" style="width: 100%; text-align: left;" 
140
|-
141
| 
142
{| style="text-align: left; margin:auto;width: 100%;" 
143
|-
144
| style="text-align: center;" | <math>\sigma _{ij}n_j -t_i^p =0 \qquad \hbox{on }\Gamma _t </math>
145
|}
146
| style="width: 5px;text-align: right;white-space: nowrap;" | (7.b)
147
|}
148
149
where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math>, <math display="inline">i=1,n_s</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.
150
151
====''Mass balance equation''====
152
153
The mass balance equation for an incompressible fluid is written as
154
155
<span id="eq-8"></span>
156
{| class="formulaSCP" style="width: 100%; text-align: left;" 
157
|-
158
| 
159
{| style="text-align: left; margin:auto;width: 100%;" 
160
|-
161
| style="text-align: center;" | <math>\varepsilon _v =0 \qquad \hbox{in }\Omega  </math>
162
|}
163
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
164
|}
165
166
Stabilized FEM techniques are needed for solving the general momentum equations in fluid mechanics when they are written in the Eulerian description. This is due to the effect of the convective acceleration terms that lead to loss of stability of standard Galerkin FEM. For Stokes flows (or for fluid flows formulated in the Lagrangian description) the convective terms vanish from the momentum equations and, consequently, the equations can be solved with the Galerkin FEM.
167
168
The problem, however, remains for obtaining stable solution for incompressible flows when an equal order interpolation is used for the velocities and the pressure, as it is the case for the element derived in the paper. This situation violates the so called inf-sup (or LBB) condition and, hence, there is a need to use  stabilization techniques for solving the mass balance equation with the FEM <span id='citeF-11'></span><span id='citeF-39'></span>[[#cite-11|[11,39]]].
169
170
==3 STABILIZED FIC FORMS OF THE MASS BALANCE EQUATION==
171
172
In our work the stabilized form of the mass balance equation is obtained using the Finite Calculus (FIC) procedure <span id='citeF-26'></span>[[#cite-26|[26]]]&#8211;<span id='citeF-31'></span>[[#cite-31|[31]]],<span id='citeF-33'></span>[[#cite-33|[33]]]&#8211;<span id='citeF-35'></span>[[#cite-35|[35]]]. We will use both the second order form of the mass balance equation for an incompressible fluid  obtained using the FIC method in space, as well as the FIC form of the mass balance equation in time. These forms have the following expressions.
173
174
===''Second order FIC mass balance equation in space''===
175
176
<span id="eq-9"></span>
177
{| class="formulaSCP" style="width: 100%; text-align: left;" 
178
|-
179
| 
180
{| style="text-align: left; margin:auto;width: 100%;" 
181
|-
182
| style="text-align: center;" | <math>\varepsilon _v + \frac{h_i^2}{12} \frac{\partial ^2 \varepsilon _v}{\partial x^2_i}=0\qquad \hbox{in }\Omega \qquad i=1,n_s </math>
183
|}
184
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
185
|}
186
187
===''FIC mass balance equation in time''===
188
189
<span id="eq-10"></span>
190
{| class="formulaSCP" style="width: 100%; text-align: left;" 
191
|-
192
| 
193
{| style="text-align: left; margin:auto;width: 100%;" 
194
|-
195
| style="text-align: center;" | <math>\varepsilon _v + \frac{\delta }{2} \frac{\partial \varepsilon _v}{\partial t}=0 \qquad \hbox{in }\Omega  </math>
196
|}
197
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
198
|}
199
200
Eq.([[#eq-9|9]]) is obtained by expressing the balance of mass in a rectangular domain of finite size (with space dimensions <math display="inline">h_1\times h_2</math> for 2D problems) and retaining higher order terms in the Taylor series expansions than those typically used in the infinitesimal theory for expressing the change of mass along the sides of the balance domain. The characteristic lengths <math display="inline">h_i</math> are related to the finite element sizes in the discretized problem, as it will be explained later.
201
202
Eq.([[#eq-10|10]]), on the other hand, is obtained by expressing the balance of mass in a domain of infinitesimal size in space and of finite dimension in time, where <math display="inline">\delta </math> is a characteristic time value.
203
204
In Eqs.([[#eq-9|9]]) and ([[#eq-10|10]]) the terms involving <math display="inline">h_i</math> and <math display="inline">\delta </math> play the role of stabilization terms respectively. The form of the <math display="inline">h_i</math> and <math display="inline">\delta </math> parameters will be defined later. Note that for <math display="inline">h_i\to 0</math> and <math display="inline">\delta \to 0</math> the standard expression of the mass balance equation ([[#eq-8|8]]) is recovered in all cases.
205
206
The derivation of Eqs.([[#eq-9|9]]) and ([[#eq-10|10]]) for a 1D problem are shown in <span id='citeF-37'></span>[[#cite-37|[37]]].
207
208
==4 FIC FORM OF THE MASS BALANCE EQUATION IN TERMS OF THE MOMENTUM EQUATIONS==
209
210
We will derive next a more useful FIC  form of the stabilized mass balance equation expressed in terms of the momentum equations.
211
212
From the momentum equations (6) we obtain (neglecting the space changes of the viscosity <math display="inline">\mu </math> in the term involving <math display="inline">\varepsilon _v</math>)
213
214
<span id="eq-11"></span>
215
{| class="formulaSCP" style="width: 100%; text-align: left;" 
216
|-
217
| 
218
{| style="text-align: left; margin:auto;width: 100%;" 
219
|-
220
| style="text-align: center;" | <math>\frac{2}{3}\mu \frac{\partial \varepsilon _v}{\partial x_i} = - \rho \frac{\partial v_i}{\partial t} + {\partial  \over \partial x_j} (2\mu \varepsilon _{ij}) + {\partial p \over \partial x_i}+b_i = -\rho {\partial v_i \over \partial t} +  \hat r_{m_i} \qquad i,j=1,n_s  </math>
221
|}
222
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
223
|}
224
225
From the last equation we deduce
226
227
<span id="eq-12"></span>
228
{| class="formulaSCP" style="width: 100%; text-align: left;" 
229
|-
230
| 
231
{| style="text-align: left; margin:auto;width: 100%;" 
232
|-
233
| style="text-align: center;" | <math>\frac{\partial \varepsilon _v}{\partial x_i} = \frac{3}{2\mu } \left[-\rho \frac{\partial v_i}{\partial t} + \hat r_{m_i}\right]\qquad i=1,n_s </math>
234
|}
235
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
236
|}
237
238
In the last two equations <math display="inline">\hat r_{m_i}</math> is a ''static momentum term'' defined as
239
240
<span id="eq-13"></span>
241
{| class="formulaSCP" style="width: 100%; text-align: left;" 
242
|-
243
| 
244
{| style="text-align: left; margin:auto;width: 100%;" 
245
|-
246
| style="text-align: center;" | <math>\hat r_{m_i}= {\partial \hat \sigma _{ij} \over \partial x_j}+b_i \quad \hbox{with}\quad  \hat \sigma _{ij}= 2\mu \varepsilon _{ij}+p \delta _{ij}   </math>
247
|}
248
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
249
|}
250
251
Let us introduce <math display="inline">\frac{\partial \varepsilon _v}{\partial x_i}</math> from Eq.([[#eq-12|12]]) intro Eq.([[#eq-9|9]]). This gives, after small algebra
252
253
<span id="eq-14"></span>
254
{| class="formulaSCP" style="width: 100%; text-align: left;" 
255
|-
256
| 
257
{| style="text-align: left; margin:auto;width: 100%;" 
258
|-
259
| style="text-align: center;" | <math>\varepsilon _v + \frac{h_i^2}{12}\frac{\partial }{\partial x_i} \left(\frac{\partial \varepsilon _v}{\partial x_i}\right)= \rho \varepsilon _v + \frac{h_i^2}{8 \mu } \frac{\partial }{\partial x_i} \left(-\rho \frac{\partial v_i}{\partial t}+\hat r_{m_i}\right)\qquad i=1,n_s </math>
260
|}
261
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
262
|}
263
264
Observing the term involving the time derivative of <math display="inline">v_i</math> in Eq.([[#eq-14|14]]) gives
265
266
<span id="eq-15"></span>
267
{| class="formulaSCP" style="width: 100%; text-align: left;" 
268
|-
269
| 
270
{| style="text-align: left; margin:auto;width: 100%;" 
271
|-
272
| style="text-align: center;" | <math>\frac{\partial }{\partial x_i}  \left(-\rho \frac{\partial v_i}{\partial t}\right)=-\rho \frac{\partial }{\partial t} \left({\partial v_i \over \partial x_i}\right)=-\rho \frac{\partial \varepsilon _v}{\partial t}  </math>
273
|}
274
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
275
|}
276
277
Substituting Eq.([[#eq-15|15]]) into ([[#eq-14|14]]) gives
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;width: 100%;" 
284
|-
285
| style="text-align: center;" | <math>\varepsilon _v + \frac{h_i^2}{8 \mu }\left(-\rho \frac{\partial \varepsilon _v}{\partial t}+ {\partial \hat r_{m_i} \over \partial x_i}   \right)=0  </math>
286
|}
287
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
288
|}
289
290
On the other hand, from Eq.([[#eq-10|10]]) we deduce
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;width: 100%;" 
297
|-
298
| style="text-align: center;" | <math>-\frac{\partial \varepsilon _v}{\partial t}=  \frac{2}{\delta } \varepsilon _v  </math>
299
|}
300
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
301
|}
302
303
Substituting Eq.([[#eq-17|17]]) into ([[#eq-16|16]]) gives
304
305
<span id="eq-18"></span>
306
{| class="formulaSCP" style="width: 100%; text-align: left;" 
307
|-
308
| 
309
{| style="text-align: left; margin:auto;width: 100%;" 
310
|-
311
| style="text-align: center;" | <math>\varepsilon _v +  \frac{h_i^2}{8\mu }\left(\frac{2\rho }{\delta } \varepsilon _v +   {\partial \hat r_{m_i} \over \partial x_i}   \right)=0  </math>
312
|}
313
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
314
|}
315
316
In the following we will assume <math display="inline">h_i = h</math> where <math display="inline">h</math> is a characteristic length that will be related to a typical average dimension of each finite element in the mesh. Multiplying Eq.([[#eq-18|18]]) by <math display="inline"> \frac{8\mu }{h^2}</math> gives, after grouping some terms,
317
318
<span id="eq-19"></span>
319
{| class="formulaSCP" style="width: 100%; text-align: left;" 
320
|-
321
| 
322
{| style="text-align: left; margin:auto;width: 100%;" 
323
|-
324
| style="text-align: center;" | <math>\displaystyle  \varepsilon _v +\tau {\partial \hat r_{m_i} \over \partial x_i} =0  </math>
325
|}
326
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
327
|}
328
329
where <math display="inline">\tau </math> is a ''stabilization parameter'' given by
330
331
<span id="eq-20"></span>
332
{| class="formulaSCP" style="width: 100%; text-align: left;" 
333
|-
334
| 
335
{| style="text-align: left; margin:auto;width: 100%;" 
336
|-
337
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{h^2}+ \frac{2\rho }{\delta } \right)^{-1}  </math>
338
|}
339
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
340
|}
341
342
Eq.([[#eq-19|19]]) is the FIC form of the stabilized mass balance equation. This equation will be taken as the starting point for deriving the stabilized FIC-FEM formulation. Note again that for <math display="inline">\tau =0</math>, the standard form of the incompressibility condition ([[#eq-8|8]]) is obtained.
343
344
==5 VARIATIONAL EQUATIONS==
345
346
===5.1 Variational expression of the momentum equation===
347
348
Multiplying Eq.(1) by arbitrary test functions <math display="inline">w_i ({x})</math> (with dimension of velocity) and integrating over the analysis domain <math display="inline">\Omega </math> gives the standard weighted residual form of the momentum equations as
349
350
<span id="eq-21"></span>
351
{| class="formulaSCP" style="width: 100%; text-align: left;" 
352
|-
353
| 
354
{| style="text-align: left; margin:auto;width: 100%;" 
355
|-
356
| style="text-align: center;" | <math>\int _\Omega w_i \left(\rho {\partial v_i \over \partial t}- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0 </math>
357
|}
358
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
359
|}
360
361
Integrating by parts  the term involving <math display="inline">\sigma _{ij}</math> and making use of the Neumann  boundary conditions ([[#eq-7.b|7.b]]), yields the standard principle of virtual power as <span id='citeF-4'></span><span id='citeF-39'></span>[[#cite-4|[4,39]]]
362
363
<span id="eq-22"></span>
364
{| class="formulaSCP" style="width: 100%; text-align: left;" 
365
|-
366
| 
367
{| style="text-align: left; margin:auto;width: 100%;" 
368
|-
369
| style="text-align: center;" | <math>\int _\Omega w_i \rho {\partial v_i \over \partial t} 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>
370
|}
371
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
372
|}
373
374
where <math display="inline">\delta \varepsilon _{ij}= \frac{1}{2}\left({\partial w_i \over \partial x_j}+{\partial w_j \over \partial x_i}\right)</math> is an arbitrary (virtual) strain rate field.
375
376
Substituting the expression of the stresses from Eq.([[#eq-2|2]]) into ([[#eq-22|22]]) gives
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;width: 100%;" 
383
|-
384
| style="text-align: center;" | <math>\int _\Omega \! w_i \rho {\partial v_i \over \partial t} d\Omega + \!\int _\Omega \left[\delta \varepsilon _{ij} 2\mu \left(\!\varepsilon _{ij} - \frac{1}{3}\varepsilon _{kk} \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>
385
|}
386
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
387
|}
388
389
Eq.([[#eq-23|23]]) can be written in matrix form as
390
391
<span id="eq-24"></span>
392
{| class="formulaSCP" style="width: 100%; text-align: left;" 
393
|-
394
| 
395
{| style="text-align: left; margin:auto;width: 100%;" 
396
|-
397
| style="text-align: center;" | <math>\displaystyle  \int _\Omega {\bf w}^T \rho {\partial {\bf v} \over \partial t} d\Omega + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf D} {\boldsymbol \varepsilon } d\Omega  + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf m} p - \int _\Omega {\bf w}^T {b}d\Omega - \!\int _{\Gamma _t} \mathbf{\bf w}^T \mathbf{\bf t}^p d\Gamma =0  </math>
398
|}
399
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
400
|}
401
402
In Eq.([[#eq-24|24]]) <math display="inline">{\bf w},{\bf v},{\boldsymbol \varepsilon }</math> are vectors containing the weighting functions, the velocities and the strain rates, respectively; <math display="inline">{\bf b}</math> and <math display="inline">\mathbf{\bf t}^p</math> are body force and surface tractions 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)
403
404
<span id="eq-25"></span>
405
{| class="formulaSCP" style="width: 100%; text-align: left;" 
406
|-
407
| 
408
{| style="text-align: left; margin:auto;width: 100%;" 
409
|-
410
| style="text-align: center;" | <math>\begin{array}{l}\mathbf{w} =[w_1,w_2,w_3]^T\quad ,\quad \mathbf{v} =[v_1,v_2,v_3]^T\\[.5cm] {\boldsymbol \varepsilon }= [\varepsilon _{11},\varepsilon _{22}, \varepsilon _{33},2\varepsilon _{12},2\varepsilon _{13}, 2\varepsilon _{23}]^T \quad ,\quad  \mathbf{\bf b} =[b_1,b_2,b_3]^T \quad ,\quad  \mathbf{\bf t}^p =[t_1^p,t_2^p,t_3^p]^T\\[.5cm]  
411
{\bf D}=\mu \left[                      \begin{array}{cccccc}\frac{4}{3} & -\frac{2}{3} &-\frac{2}{3} & 0 & 0 & 0 \\                         &   \frac{4}{3} & -\frac{2}{3} &0  &0  & 0 \\                         &  &  \frac{4}{3} & 0 & 0 & 0 \\                         &  &  & 1 & 0 & 0 \\                        {Sym.} &  &  &  & 1 & 0 \\                         &  &  &  &  & 1 \\                      \end{array} \right]\qquad ,\qquad {\bf m}=[1,1,1,0,0,0]^T \end{array} </math>
412
|}
413
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
414
|}
415
416
The 2D form of above expressions is straightforward. For instance, the 2D expression of matrix <math display="inline"> {\bf D}</math> is obtained by deleting the rows and columns 3, 5 and 6 in Eq.([[#eq-25|25]]).
417
418
<br/>
419
420
'''Remark 1.''' From the definition of <math display="inline">{\bf m}</math> and <math display="inline">{\bf D}</math> and Eqs.([[#eq-2|2]]), ([[#eq-3|3]]) and ([[#eq-5|5]])  we deduce
421
422
{| class="formulaSCP" style="width: 100%; text-align: left;" 
423
|-
424
| 
425
{| style="text-align: left; margin:auto;width: 100%;" 
426
|-
427
| style="text-align: center;" | <math>{\boldsymbol \sigma } =  {\bf s}  + {\bf m}p   \quad ,\quad   {\bf s}= {\bf D}{\boldsymbol \varepsilon }  \quad \hbox{and}\quad \varepsilon _v= {\bf m}^T {\boldsymbol \varepsilon } </math>
428
|}
429
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
430
|}
431
432
where <math display="inline">{\boldsymbol \sigma } =[\sigma _{11},\sigma _{22},\sigma _{33}, \sigma _{12},\sigma _{13}, \sigma _{31}]^T,~{s}=[s_{11},s_{22},s_{33}, s_{12},s_{13}, s_{23}]^T</math> are the stress and deviatoric stress vectors, respectively.<br/>
433
434
===5.2 Variational expression of the stabilized mass balance equation===
435
436
Multiplying equation ([[#eq-19|19]]) by arbitrary test functions <math display="inline">q(\mathbf{x})</math> (with dimension of pressure) defined over the analysis domain <math display="inline">\Omega </math> and integrating over <math display="inline">\Omega </math> gives
437
438
<span id="eq-27"></span>
439
{| class="formulaSCP" style="width: 100%; text-align: left;" 
440
|-
441
| 
442
{| style="text-align: left; margin:auto;width: 100%;" 
443
|-
444
| style="text-align: center;" | <math>\int _\Omega q \varepsilon _v d\Omega  + \int _\Omega q \tau {\partial \hat r_{m_i} \over \partial x_i} d\Omega =0  </math>
445
|}
446
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
447
|}
448
449
Let us integrate by parts the second integral in Eq.([[#eq-27|27]]) over a mesh of <math display="inline">N_e</math> elements each one  with area <math display="inline">\Omega ^e </math> and boundary <math display="inline">\Gamma ^e</math>. This yields
450
451
<span id="eq-28"></span>
452
{| class="formulaSCP" style="width: 100%; text-align: left;" 
453
|-
454
| 
455
{| style="text-align: left; margin:auto;width: 100%;" 
456
|-
457
| style="text-align: center;" | <math>\sum \limits _{e=1}^{N_e} \left\{\int _{\Omega ^e} q \varepsilon _v d\Omega - \int _{\Omega ^e} \tau {\partial q \over \partial x_i} \hat{r}_{m_{i}} d \Omega  +\int _{\Gamma ^e} q \tau  n_i \hat{r}_{m_{i}} d \Gamma \right\}=0   </math>
458
|}
459
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
460
|}
461
462
where <math display="inline">n_i</math> are the components of the unit normal vector to the element boundary <math display="inline">\Gamma ^e</math>. Note that <math display="inline">\Gamma ^e</math> is a side for a 3-noded triangle and a face for a 4-noded tetrahedron.
463
464
<br/>
465
466
'''Remark 2.''' In Eq.([[#eq-28|28]]), the  second integral over the domain <math display="inline">\Omega ^e</math> involving <math display="inline">{\partial q \over \partial x_i}</math> is zero for a piecewise constant approximation of the weighting function <math display="inline">q</math> over the  elements. This situation occurs for the P1/P0+ elements presented in this work.
467
468
From the definition of <math display="inline">\hat{r}_{m_{i}}</math> of Eqs.(13) and (1) we deduce
469
470
<span id="eq-29"></span>
471
{| class="formulaSCP" style="width: 100%; text-align: left;" 
472
|-
473
| 
474
{| style="text-align: left; margin:auto;width: 100%;" 
475
|-
476
| style="text-align: center;" | <math>\hat{r}_{m_{i}} = \rho {\partial v_i \over \partial t}+ \frac{2\mu }{3}{\partial \varepsilon _v \over \partial x_i} </math>
477
|}
478
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
479
|}
480
481
Multiplying by the normal components gives
482
483
<span id="eq-30"></span>
484
{| class="formulaSCP" style="width: 100%; text-align: left;" 
485
|-
486
| 
487
{| style="text-align: left; margin:auto;width: 100%;" 
488
|-
489
| style="text-align: center;" | <math>\hat{r}_{m_{i}}n_i = \rho {\partial v_n \over \partial t}+\frac{2\mu }{3}{\partial \varepsilon _v \over \partial n} </math>
490
|}
491
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
492
|}
493
494
where <math display="inline">v_n =v_in_i</math> is the projection of the velocity in the direction of the normal to the element boundary and <math display="inline">{\partial \varepsilon _v \over \partial n}={\partial \varepsilon _v \over \partial x_i}n_i</math> is the derivative of <math display="inline">\varepsilon _v</math> in the normal direction.
495
496
The term <math display="inline">\mu {\partial \varepsilon _v \over \partial n}</math> at an element boundary <math display="inline">\Gamma _{ea}</math> connecting elements <math display="inline">e</math> and <math display="inline">a</math> can be computed as
497
498
<span id="eq-31"></span>
499
{| class="formulaSCP" style="width: 100%; text-align: left;" 
500
|-
501
| 
502
{| style="text-align: left; margin:auto;width: 100%;" 
503
|-
504
| style="text-align: center;" | <math>\frac{2}{3}\mu {\partial \varepsilon _v \over \partial n} = \frac{2}{l_{ea}} \left(\frac{2}{3}\mu ^+ \varepsilon _v^+ - \frac{2}{3}\mu ^- \varepsilon _v^-   \right)  </math>
505
|}
506
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
507
|}
508
509
where <math display="inline">(\cdot )^+</math> and <math display="inline">(\cdot )^-</math> denote the values of the relevant magnitude at the external and internal points adjacent to the element boundary, respectively and <math display="inline">l_{ea}</math> is the characteristic length of the boundary. For triangles <math display="inline">l_{ea}</math> is taken as the side lenght  (Figure [[#img-1|1]]). For 4-noded tetrahedra <math display="inline">l_{ea} = 2 (\Omega _{ea})^{1/2}</math> where <math display="inline">\Omega _{ea}</math> is the area of the triangular face connecting elements <math display="inline">e</math> and <math display="inline">a</math>.
510
511
<div id='img-1a'></div>
512
<div id='img-1b'></div>
513
<div id='img-1'></div>
514
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
515
|-style="text-align: center; font-size: 75%;"
516
| (a) 
517
| (b)
518
|
519
|-
520
|[[Image:Draft_Samper_536119838-Fig1a.png|300px|]]
521
|[[Image:Draft_Samper_536119838-Fig1b.png|300px|]]
522
| [[File:Draft_Samper_536119838_2271_formulafig1.jpg]]
523
|- style="text-align: center; font-size: 75%;"
524
| colspan="2" |(c)
525
|
526
|-
527
|colspan="2" | [[Image:Draft_Samper_536119838-Fig1c.png|400px|]]
528
|
529
|- style="text-align: center; font-size: 75%;"
530
| colspan="3" | '''Figure 1:''' (a) Patch of four triangles associated to a central 3-noded triangle <math>e</math>. (b) Definition of <math>{\partial \varepsilon _v \over \partial n} </math> at the side <math>ij</math> of element <math>e</math>. (c) Equilibrium of tractions at the side <math>ij</math> adjacent to elements <math>e</math> and <math>a</math>. (d) Jump of <math>\frac{2}{3}\mu \varepsilon _v</math> across a side of element <math>e</math>
531
|}
532
533
Eq.([[#eq-31|31]]) is just a simple (yet effective) model for computing <math display="inline">\mu {\partial \varepsilon _v \over \partial n}</math> at the element boundary. Indeed, more sophisticated procedures can be used <span id='citeF-9'></span>[[#cite-9|[9]]].
534
535
From the Neumann boundary conditions (7b) and Eq.([[#eq-31|31]]) we deduce
536
537
{| class="formulaSCP" style="width: 100%; text-align: left;" 
538
|-
539
| 
540
{| style="text-align: left; margin:auto;width: 100%;" 
541
|-
542
| style="text-align: center;" | <math>\frac{2}{3}\mu ^+ \varepsilon _v^+  = 2\mu ^+  {\partial v_n^+ \over \partial n} + p^+ - t^+_n \quad \hbox{ at }\Gamma _{ea}^+</math>
543
| style="width: 5px;text-align: right;white-space: nowrap;" | (32.a)
544
|-
545
| style="text-align: center;" | <math>  \frac{2}{3}\mu ^- \varepsilon _v^-  = 2 \mu ^- {\partial v_n^- \over \partial n} + p^- - t^-_n \quad \hbox{ at }\Gamma _{ea}^- </math>
546
| style="width: 5px;text-align: right;white-space: nowrap;" | (32.b)
547
|}
548
|}
549
550
where indices <math display="inline">+</math> and &#8211; denote values at the external and internal sides <math display="inline">\Gamma _{ea}^+</math> and <math display="inline">\Gamma _{ea}^-</math> of the element boundary, respectively. In Eqs.([[#5.2 Variational expression of the stabilized mass balance equation|5.2]]) <math display="inline">t_n^+</math> and <math display="inline">t_n^-</math> are the normal tractions respectively acting at <math display="inline">\Gamma _{ea}^+</math> and <math display="inline">\Gamma _{ea}^-</math> (Figure [[#img-1|1]]c).
551
552
Substituting Eqs.([[#5.2 Variational expression of the stabilized mass balance equation|5.2]]) into ([[#eq-31|31]]) gives
553
554
<span id="eq-33"></span>
555
{| class="formulaSCP" style="width: 100%; text-align: left;" 
556
|-
557
| 
558
{| style="text-align: left; margin:auto;width: 100%;" 
559
|-
560
| style="text-align: center;" | <math>\frac{2}{3}\mu {\partial \varepsilon _v \over \partial n} = \frac{2}{l_{ea}} \left[\!\!\left[2\mu {\partial v_n \over \partial n}+p  -t_n \right]\!\!\right]_{\Gamma ^e}  </math>
561
|}
562
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
563
|}
564
565
where <math display="inline">\left[\!\left[a   \right]\!\right]</math> denotes the jump of the magnitude <math display="inline">a</math> across the element boundary <math display="inline">\Gamma ^e</math>, i.e. <math display="inline">\left[\!\left[a   \right]\!\right]= a^+-a^-</math>.
566
567
From the equilibrium of tractions at the element boundaries we obtain
568
569
<span id="eq-34"></span>
570
{| class="formulaSCP" style="width: 100%; text-align: left;" 
571
|-
572
| 
573
{| style="text-align: left; margin:auto;width: 100%;" 
574
|-
575
| style="text-align: center;" | <math>\left[\!\left[t_n \right]\!\right]= t_n^+ - t_n^- = t_n^p = \left[{t}^p \right]^T {n}^+  </math>
576
|}
577
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
578
|}
579
580
where <math display="inline">t_n^p</math> is the external normal traction acting on the element side in the direction of the  normal to <math display="inline">\Gamma ^+_{ea}</math> (Figure [[#img-1|1]]). Clearly ''for unloaded boundaries'' <math display="inline">t_n^p=0</math> and, consequently, <math display="inline">\left[\!\left[t_n \right]\!\right]=0</math>.
581
582
We note that above derivations are independent of the choice of <math display="inline">\Gamma _{ea}^+</math> and <math display="inline">\Gamma _{ea}^-</math>.
583
584
Substituting Eq.([[#eq-33|33]]) into ([[#eq-30|30]]) and this one into ([[#eq-28|28]]) gives the final stabilized variational form of the mass balance equation as
585
586
<span id="eq-35"></span>
587
{| class="formulaSCP" style="width: 100%; text-align: left;" 
588
|-
589
| 
590
{| style="text-align: left; margin:auto;width: 100%;" 
591
|-
592
| style="text-align: center;" | <math>\sum \limits _{e=1}^{N_e} \left\{\int _{\Omega ^e} q \varepsilon _v d\Omega - \int _{\Omega ^e} \tau  {\partial q \over \partial x_i} \hat r_{m_i} d\Omega + \int _{\Gamma ^e} \frac{2\tau }{l^e}q \left(\frac{\rho l^e}{2} {\partial v_n \over \partial t} +\left[\!\!\left[2\mu {\partial v_n \over \partial n}+p -t_n  \right]\!\!\right]\right)d\Gamma \right\}=0  </math>
593
|}
594
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
595
|}
596
597
where <math display="inline">l^e</math> are the lengths of the element sides.
598
599
For element boundaries laying on the Dirichlet boundary <math display="inline">\Gamma _v</math>, <math display="inline">\hat r_{m_i} n_i =0</math> and the boundary term vanishes from Eq.([[#eq-35|35]]). This is justified by the fact that both <math display="inline">{\partial v_n \over \partial t}</math> and <math display="inline">{\partial \varepsilon _v \over \partial n}</math> are zero at <math display="inline">\Gamma _v</math>.
600
601
'''Remark 3.''' At element boundaries laying on an external Neumann boundary <math display="inline">\Gamma _t</math>, <math display="inline">t^+_n=0</math>, <math display="inline">{n}^+=-{n}^- = {n}^e</math> and, hence, <math display="inline">\left[\!\left[t_n \right]\!\right]=-t_n^p</math> at <math display="inline">\Gamma _t</math>.
602
603
'''Remark 4.''' The boundary integral in Eq.([[#eq-35|35]]) resembles the jump stabilization terms across element sides proposed by different authors for fluid flow problems <span id='citeF-3'></span><span id='citeF-9'></span><span id='citeF-12'></span><span id='citeF-14'></span><span id='citeF-17'></span><span id='citeF-23'></span>[[#cite-3|[3,9,12,14,17,23]]]. It is remarkable that Eq.([[#eq-35|35]]) emanates naturally from the FIC formulation. Also, the form of Eq.([[#eq-35|35]]) introduces the effect of the acceleration in the normal direction to the element side, as well as the effect of the external surface tractions acting on the element side. These two terms have proven to be important for the improved mass conservation and overall stability of the numerical solution of the multifluid problems solved in this paper and also for free surface homogeneous fluid flows <span id='citeF-37'></span>[[#cite-37|[37]]].
604
605
'''Remark 5.''' For steady state problems the acceleration term <math display="inline">\rho {\partial v_n \over \partial t}</math> vanishes from the boundary integral in Eq.([[#eq-35|35]]).
606
607
==6 FEM DISCRETIZATION==
608
609
We discretize  the analysis domain into a mesh of <math display="inline">N_e</math> finite elements in the standard manner. In our work we will choose simple 3-noded triangles (for 2D problems) and 4-noded tetrahedra (for 3D problems). The velocity is linearly interpolated in terms of the nodal values, while a discontinuous linear interpolation for the pressure over each element is chosen, i.e.
610
611
<span id="eq-36"></span>
612
{| class="formulaSCP" style="width: 100%; text-align: left;" 
613
|-
614
| 
615
{| style="text-align: left; margin:auto;width: 100%;" 
616
|-
617
| style="text-align: center;" | <math>v_i = \sum \limits _{j=1}^n N_j \bar v_i^j \quad ,\quad p = \bar p^e + \tilde p({\bf x})      </math>
618
|}
619
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
620
|}
621
622
where <math display="inline">N_i</math> are the standard linear shape functions for simplicial elements, <math display="inline">n</math> is the number of element nodes (<math display="inline">n=3/4</math> for 2D/3D problems), <math display="inline">\bar{v}_i^j</math> denotes the value of the <math display="inline">i</math>th velocity component for the <math display="inline">j</math>th node of an element, <math display="inline">\bar p^e</math> is a constant pressure field over each element and <math display="inline"> \tilde {p}({\bf x})</math> is a discontinuous linear pressure field chosen as (Figure [[#img-2|2]])
623
624
<span id="eq-37"></span>
625
{| class="formulaSCP" style="width: 100%; text-align: left;" 
626
|-
627
| 
628
{| style="text-align: left; margin:auto;width: 100%;" 
629
|-
630
| style="text-align: center;" | <math>\tilde {p}({\bf x}) = ({\bf x}_{c} - {\bf x})^{T} {\bf b} </math>
631
|}
632
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
633
|}
634
635
where <math display="inline">{x}_{c}</math> are the coordinates of the element midpoint.
636
637
<div id='img-2'></div>
638
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
639
|-
640
|[[Image:Draft_Samper_536119838-discontinuous-pressure.png|540px|One dimensional representation of the discontinuous pressure field. (a) Constant element pressure (̄p<sup>e</sup>) and linear pressure field (̃p). (b) Pressure field p=̄p<sup>e</sup>+ ̃p over the element. (c) Discontinuous linear pressure field over a three 1D element patch]]
641
|- style="text-align: center; font-size: 75%;"
642
| colspan="1" | '''Figure 2:''' One dimensional representation of the discontinuous pressure field. (a) Constant element pressure (<math>\bar p^e</math>) and linear pressure field (<math>\tilde p</math>). (b) Pressure field <math>p=\bar p^e + \tilde p</math> over the element. (c) Discontinuous linear pressure field over a three 1D element patch
643
|}
644
645
Note that <math display="inline">\tilde {p}({x})</math> satisfies the steady state form of the momentum equation for a linear velocity field and  constant body forces <math display="inline">b_i</math>, i.e.
646
647
<span id="eq-38"></span>
648
{| class="formulaSCP" style="width: 100%; text-align: left;" 
649
|-
650
| 
651
{| style="text-align: left; margin:auto;width: 100%;" 
652
|-
653
| style="text-align: center;" | <math>{\partial \tilde p \over \partial x_i}+b_i=0 </math>
654
|}
655
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
656
|}
657
658
The actual unknowns of the problem are the nodal velocities <math display="inline">\bar v_i^j</math> and the constant pressure field <math display="inline">\bar p^e</math> for each element (Figure [[#img-3|3]]).
659
660
We will choose now <math display="inline">N_e</math> piecewise constant unit weighting functions <math display="inline">q^e ({\bf x})</math>  so that <math display="inline">q^e ({\bf x})=1</math> and <math display="inline"> {\partial q^e \over \partial x_i}=0 </math> if <math display="inline"> {x} \in \Omega ^e</math> and <math display="inline">q^e ({\bf x}) =0</math> if <math display="inline"> {x} \not \in \Omega ^e</math>. With these assumptions the variational form ([[#eq-35|35]]) for the stabilized mass balance equation simplifies for 3-noded triangles to
661
662
<span id="eq-39"></span>
663
{| class="formulaSCP" style="width: 100%; text-align: left;" 
664
|-
665
| 
666
{| style="text-align: left; margin:auto;width: 100%;" 
667
|-
668
| style="text-align: center;" | <math>\int _{\Omega ^e}  \varepsilon _v d\Omega + \sum \limits _{i=1}^{3}  2 \tau \left( \frac{\rho l_i^e}{2} {\partial v_n \over \partial t}+ \left[\!\!\left[2\mu  {\partial v_n \over \partial n}+p -t_n\right]\!\!\right] \right) =0 \quad ,\quad e=1,N_e  </math>
669
|}
670
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
671
|}
672
673
where the sum in the second term extends over the three sides of each triangular element and <math display="inline">l_i^e</math> is the length of the <math display="inline">i</math>th element side.
674
675
<div id='img-3'></div>
676
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
677
|-
678
|[[Image:Draft_Samper_536119838-nodal-velocities.png|540px|Nodal velocities and element pressure variables at a patch of four 3-noded triangles]]
679
|- style="text-align: center; font-size: 75%;"
680
| colspan="1" | '''Figure 3:''' Nodal velocities and element pressure variables at a patch of four 3-noded triangles
681
|}
682
683
'''Remark 6.''' For 4-noded tetrahedra the sum in Eq.([[#eq-39|39]]) extends over the four faces of the tetrahedra  and <math display="inline">l_i^e</math> is a characteristic distance for the <math display="inline">i</math>th face with area <math display="inline">\Omega _i^e</math> computed as <math display="inline">l_i^e = 2(\Omega _i^e)^{1/2}</math> where <math display="inline">\Omega _i^e</math> is the area of the face.
684
685
Substituting the approximation ([[#eq-36|36]]) into Eqs.([[#eq-39|39]]) and ([[#eq-24|24]]) and choosing a Galerkin form with <math display="inline">w_j =N_j</math> gives the discretized expression of the momentum and (stabilized) mass balance equations as
686
687
<span id="eq-40.a"></span>
688
{| class="formulaSCP" style="width: 100%; text-align: left;" 
689
|-
690
| 
691
{| style="text-align: left; margin:auto;width: 100%;" 
692
|-
693
| style="text-align: center;" | <math>{\bf M} {\dot{\bar{\bf v}}} + {\bf K}\bar{\bf v}+{\bf Q}\bar {\bf p}- {\bf f}_v={\bf 0} </math>
694
|}
695
| style="width: 5px;text-align: right;white-space: nowrap;" | (40.a)
696
|}
697
698
<span id="eq-40.b"></span>
699
{| class="formulaSCP" style="width: 100%; text-align: left;" 
700
|-
701
| 
702
{| style="text-align: left; margin:auto;width: 100%;" 
703
|-
704
| style="text-align: center;" | <math>{\bf Q}^T \bar{\bf v} + {\bf S}\bar{\bf p}-{\bf f}_p={\bf 0} </math>
705
|}
706
| style="width: 5px;text-align: right;white-space: nowrap;" | (40.b)
707
|}
708
709
with
710
711
<span id="eq-40.c"></span>
712
{| class="formulaSCP" style="width: 100%; text-align: left;" 
713
|-
714
| 
715
{| style="text-align: left; margin:auto;width: 100%;" 
716
|-
717
| style="text-align: center;" | <math>\bar{\bf v} = \left\{\begin{matrix}\bar{\bf v}^1\\\bar{\bf v}^2\\\vdots \\ \bar{\bf v}^N\end{matrix}\right\}~~,~~ \bar{\bf v}^j = [\bar{\bf v}^j_1,\bar{\bf v}^j_2,\bar{\bf v}^j_3]^T \quad \hbox{and } \bar{\bf p} = [\bar{p}^1,\bar{p}^2,\cdots , \bar{p}^{N_e}]^{T} </math>
718
|}
719
| style="width: 5px;text-align: right;white-space: nowrap;" | (40.c)
720
|}
721
722
where <math display="inline">N</math> is the number of nodes in the mesh. The different matrices and vectors in Eqs.([[#6 FEM DISCRETIZATION|6]]a,b) are assembled from the element contributions given in Box 1 for 3-noded triangles.
723
724
<div class="center">
725
{| class="formulaSCP" style="width: 80%; text-align: center;" 
726
|-style="text-align: center;
727
| 
728
{| style="text-align: center; margin:auto;width: 80%;font-size:85%;" 
729
|<math>\begin{array}{l}
730
\displaystyle \textbf{M}_{ij}^e  =\int_{\Omega^e}\rho \textbf{N}_i^T \textbf{N}_jd\Omega ~,~\textbf{K}^e_{ij} =\int_{\Omega^e} \textbf{B}_i^T \textbf{D}\textbf{B}_j d\Omega ~,~\textbf{Q}^e_{ij}  =\int_{\Omega^e} \textbf{B}_i^T \textbf{m}d\Omega\\[.4cm]
731
\displaystyle {S}_{ee} = - 2 \left(\tau^{ea} + \tau^{eb}+\tau^{ec}\right)\quad ,\quad {S}_{eb} = 2\tau^{eb}~,~{S}_{ec} = 2\tau^{ec}~,~{S}_{ea} = 2\tau^{ea}\\[.4cm]
732
\displaystyle{\bf B}_i = \left[ \begin{matrix}
733
\displaystyle\frac{N_i}{x_1} &0\\
734
0&  \displaystyle\frac{N_i}{x_2}\\
735
\displaystyle\frac{N_i}{x_2}&\displaystyle\frac{N_i}{x_1}\\[.25cm]
736
\end{matrix}  \right]\quad , \quad
737
{\bf N}_i = \left[
738
              \begin{array}{ccc}
739
                N_i & 0  \\
740
                0 & N_i  \\
741
              \end{array}
742
            \right] \quad \quad , \quad {\bf m} = [1,1,0]^T\\\\
743
\displaystyle \textbf{f}^e_{v_i}=  \int_{\Omega^e}\textbf{N}_i^T \textbf{b} d\Omega + \int_{\Gamma_t} \textbf{N}_i^T {\bf t}d\Gamma\\[.4cm]
744
\displaystyle f_{p_e}=-\sum\limits_{r=a,b,c} 2\tau^{er} \left( \frac{\rho l_{er}}{2} \frac{v_n}{t}+
745
\left[\!\!\left[ 2\mu \frac{v_n}{n}+ ({\bf x}_c^T - {\bf x}_{er}^T)  {\bf b}-{t}_n  \right]\!\!\right]\right)
746
\end{array}</math>
747
|}
748
{|
749
|-style="text-align: left; margin:auto;width: 80%;font-size:85%;" 
750
|where <math>l_{er}</math> is the length of the side connecting elements ''e'' and ''r'' and <math>{\bf x}_{er}</math>  is the mid-point of the side. For continuous body forces between elements <math>\left[\!\left[{\bf x}_{er}^T {\bf b} \right]\!\right]=0</math>
751
|-style="text-align: left; margin:auto;width: 80%;font-size:85%;" 
752
|The expression of <math>\tau^{er}</math> is given in Eqs.(45) and (46).
753
|}
754
|}
755
</div>
756
757
<div class="center" style="font-size: 75%;">
758
'''Box 1'''. Element matrices and vectors in Eqs.(\ref{eq36})  for 3-noded triangles (Figure 1)</div>
759
760
761
762
'''Remark 7.''' The FEM approximation ([[#eq-36|36]]) yields the strain rates and the stresses within an element in terms of the nodal velocities and pressure as
763
764
<span id="eq-41"></span>
765
{| class="formulaSCP" style="width: 100%; text-align: left;" 
766
|-
767
| 
768
{| style="text-align: left; margin:auto;width: 100%;" 
769
|-
770
| style="text-align: center;" | <math>\begin{array}{c}{\boldsymbol \varepsilon } = \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j \quad ,\quad \varepsilon _v = \mathbf{m}^T \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j\\ \displaystyle \mathbf{s} = \mathbf{D} \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j \quad ,\quad  {\boldsymbol \sigma } =  \mathbf{D} \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j + \mathbf{m} \bar p^e\end{array}  </math>
771
|}
772
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
773
|}
774
775
where <math display="inline">\mathbf{B}_j</math> is given in Box 1 for 2D problems.
776
777
==7 SOLUTION OF THE DISCRETIZED EQUATIONS==
778
779
The discretized form in time of the system of Eqs.([[#6 FEM DISCRETIZATION|6]]) is expressed as
780
781
<span id="eq-42.a"></span>
782
{| class="formulaSCP" style="width: 100%; text-align: left;" 
783
|-
784
| 
785
{| style="text-align: left; margin:auto;width: 100%;" 
786
|-
787
| style="text-align: center;" | <math>{\bf M} \frac{{}^{n+1}\bar{\bf v}-{}^{n}\bar{\bf v}}{\Delta t} + {\bf K} {}^{n+1}\bar{\bf v}+{\bf Q}{}^{n+1}\bar {\bf p}-{\bf f}_v={\bf 0} </math>
788
|}
789
| style="width: 5px;text-align: right;white-space: nowrap;" | (42.a)
790
|}
791
792
<span id="eq-42.b"></span>
793
{| class="formulaSCP" style="width: 100%; text-align: left;" 
794
|-
795
| 
796
{| style="text-align: left; margin:auto;width: 100%;" 
797
|-
798
| style="text-align: center;" | <math>{\bf Q}^T  {}^{n+1}\bar{\bf v}+ {\bf S} {}^{n+1}\bar{\bf p}-{\bf f}_p={\bf 0} </math>
799
|}
800
| style="width: 5px;text-align: right;white-space: nowrap;" | (42.b)
801
|}
802
803
A symmetric monolithic form of Eq.([[#7 SOLUTION OF THE DISCRETIZED EQUATIONS|7]]) can be written as
804
805
<span id="eq-43"></span>
806
{| class="formulaSCP" style="width: 100%; text-align: left;" 
807
|-
808
| 
809
{| style="text-align: left; margin:auto;width: 100%;" 
810
|-
811
| style="text-align: center;" | <math>\begin{bmatrix}\frac{1}{\Delta t}{\bf M}+ {\bf K}& {\bf Q} \\[.2cm]  {\bf Q}^T  &  {\bf S} \end{bmatrix} \left\{\begin{matrix}{}^{n+1}\bar{\bf v} \\[.2cm]{}^{n+1}\bar {\bf p}\\[.2cm]\end{matrix}\right\}=\left\{\begin{matrix}{\bf f}_v+   \frac{1}{\Delta t} {\bf M}  {}^{n}\bar{\bf v}\\[.2cm]  {\bf f}_p  \end{matrix}\right\} </math>
812
|}
813
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
814
|}
815
816
The steady state form of Eq.([[#eq-43|43]]) is simply
817
818
<span id="eq-44"></span>
819
{| class="formulaSCP" style="width: 100%; text-align: left;" 
820
|-
821
| 
822
{| style="text-align: left; margin:auto;width: 100%;" 
823
|-
824
| style="text-align: center;" | <math>\begin{bmatrix}{\bf K}&{\bf Q}\\[.2cm]{\bf Q}^T  &  {\bf S} \end{bmatrix} \left\{\begin{matrix}\bar{v} \\[.2cm]\bar {\bf p} \end{matrix}\right\}= \left\{\begin{matrix}{\bf f}_v\\[.2cm]{\bf f}_p  \end{matrix}\right\} </math>
825
|}
826
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
827
|}
828
829
Note that the terms involving the normal velocity to the element, the normal tractions to the side and the term emanating form the side discontinuous  linear pressure field <math display="inline">\tilde p</math> have been incorporated into the force vector <math display="inline">{\bf f}_p</math> in Eq.([[#eq-40.b|40.b]]). An obvious alternative will be to add the contribution of these terms to matrices <math display="inline">\bf Q</math> and <math display="inline">\bf S</math> in the l.h.s. of this equation.
830
831
'''Remark 8.''' The stabilization parameter <math display="inline">\tau </math> of Eq.([[#eq-20|20]]) is computed at each element boundary <math display="inline">ij</math> connecting elements <math display="inline">i</math> and <math display="inline">j</math> as
832
833
<span id="eq-45.a"></span>
834
{| class="formulaSCP" style="width: 100%; text-align: left;" 
835
|-
836
| 
837
{| style="text-align: left; margin:auto;width: 100%;" 
838
|-
839
| style="text-align: center;" | <math>\tau = \tau ^{ij}= \left(\frac{8\mu _{ij}}{l^2_{ij}}+\frac{2\rho _{ij}}{\Delta t}   \right)^{-1} </math>
840
|}
841
| style="width: 5px;text-align: right;white-space: nowrap;" | (45.a)
842
|}
843
844
where <math display="inline">l_{ij}</math> is a characteristic length of the element boundary. For triangles <math display="inline">l_{ij}</math> is taken as the side length. For 4-noded tetrahedra <math display="inline">l_{ij}=2 (\Omega _{ij})^{1/2}</math>, where <math display="inline">\Omega _{ij}</math> is the area of the face connecting elements <math display="inline">i</math> and <math display="inline">j</math>. The material parameters <math display="inline">\mu _{ij}</math> and <math display="inline">\rho _{ij}</math> are computed as
845
846
<span id="eq-45.b"></span>
847
{| class="formulaSCP" style="width: 100%; text-align: left;" 
848
|-
849
| 
850
{| style="text-align: left; margin:auto;width: 100%;" 
851
|-
852
| style="text-align: center;" | <math>\mu _{ij} = \frac{1}{2} (\mu _i + \mu _j) \quad ,\quad \rho _{ij}=  \frac{1}{2} (\rho _i + \rho _j) </math>
853
|}
854
| style="width: 5px;text-align: right;white-space: nowrap;" | (45.b)
855
|}
856
857
where indices <math display="inline">i</math> and <math display="inline">j</math> denote the element where the material parameter is computed.
858
859
'''Remark 9.''' For ''the steady state case'' the stabilization parameter of Eq.([[#eq-45.a|45.a]]) is computed as
860
861
<span id="eq-46"></span>
862
{| class="formulaSCP" style="width: 100%; text-align: left;" 
863
|-
864
| 
865
{| style="text-align: left; margin:auto;width: 100%;" 
866
|-
867
| style="text-align: center;" | <math>\tau ^{ij}= \left(\frac{8\mu _{ij}}{l^2_{ij}}+\frac{2\rho _{ij} \vert{v}_{ij}\vert }{l_{ij}}   \right)^{-1} </math>
868
|}
869
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
870
|}
871
872
with <math display="inline"> {v}_{ij}</math> being the velocity vector of the mid-point of the boundary connecting elements <math display="inline">i</math> and <math display="inline">j</math>.
873
874
The above expressions of the stabilization parameter are similar to those used by many stabilization methods <span id='citeF-11'></span><span id='citeF-39'></span>[[#cite-11|[11,39]]].
875
876
==8 EXTENSION TO THE NAVIER-STOKES EQUATIONS==
877
878
The momentum equations for Navier-Stokes flows are written in the Eulerian framework as
879
880
<span id="eq-47"></span>
881
{| class="formulaSCP" style="width: 100%; text-align: left;" 
882
|-
883
| 
884
{| style="text-align: left; margin:auto;width: 100%;" 
885
|-
886
| style="text-align: center;" | <math>r_{m_i}:= \rho \left(\frac{\partial v_i}{\partial t}+v_j{\partial v_{i} \over \partial x_j} \right)- {\partial \sigma _{ij} \over \partial x_j}-b_i=0  </math>
887
|}
888
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
889
|}
890
891
Note that the difference with Eq.([[#eq-1|1]]) is that the convective acceleration terms are now accounted for in Eq.([[#eq-47|47]]).
892
893
The constitutive equation and the boundary conditions are identical to Eqs.([[#eq-2|2]])&#8211;([[#eq-4|4]]) and (7), respectively.
894
895
The mass balance equation is given by Eq.([[#eq-8|8]]). The stabilized form is given by Eq.([[#eq-19|19]]).
896
897
The FEM  solution of the Navier-Stokes equation for incompressible fluids requires a stabilized numerical method that can capture the internal layers introduced by the convective acceleration terms, as well as for satisfying the inf-sup condition introduced by the incompressibility constraint. Different stabilization procedures have been proposed in the past two decades. An overview can be found in <span id='citeF-11'></span><span id='citeF-39'></span>[[#cite-11|[11,39]]]. In our work we use the FIC approach for obtaining stabilized FEM solutions for the Navier-Stokes equations <span id='citeF-26'></span>[[#cite-26|[26]]]&#8211;<span id='citeF-31'></span>[[#cite-31|[31]]],<span id='citeF-33'></span>[[#cite-33|[33]]]&#8211;<span id='citeF-35'></span>[[#cite-35|[35]]].
898
899
In the FIC method the momentum equations are derived in a balance domain of finite size. This yields a non-local form of the equations. For the <math display="inline">i</math>th momentum equation we have
900
901
<span id="eq-48"></span>
902
{| class="formulaSCP" style="width: 100%; text-align: left;" 
903
|-
904
| 
905
{| style="text-align: left; margin:auto;width: 100%;" 
906
|-
907
| style="text-align: center;" | <math>r_{m_i} - \frac{h_{ij}}{2} \frac{\partial r_{m_i}}{\partial x_j}=0 \quad i,j=1,n_s \hbox{ (sum in }j\hbox{ only)} </math>
908
|}
909
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
910
|}
911
912
where <math display="inline">h_{ij}</math> are characteristic distances that define the balance domain <span id='citeF-28'></span>[[#cite-28|[28]]]&#8211;<span id='citeF-30'></span>[[#cite-30|[30]]].
913
914
For consistency reasons, in the FIC method  the Neumann boundary conditions expressing balance of tractions at the boundary <math display="inline">\Gamma _t</math> are also derived in a finite domain adjacent to the boundary. The modified Neumann boundary conditions are written as
915
916
<span id="eq-49"></span>
917
{| class="formulaSCP" style="width: 100%; text-align: left;" 
918
|-
919
| 
920
{| style="text-align: left; margin:auto;width: 100%;" 
921
|-
922
| style="text-align: center;" | <math>\sigma _{ij}n_h - t_i^p + \frac{h_{in}}{2} r_{m_i} =0 \quad i,j=1,n_s \quad \hbox{at } \Gamma _t  </math>
923
|}
924
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
925
|}
926
927
where <math display="inline">h_{in} = h_{ij} n_j</math>.
928
929
Applying the standard weighted residual procedure to Eq.([[#eq-48|48]]) we obtain
930
931
<span id="eq-50"></span>
932
{| class="formulaSCP" style="width: 100%; text-align: left;" 
933
|-
934
| 
935
{| style="text-align: left; margin:auto;width: 100%;" 
936
|-
937
| style="text-align: center;" | <math>\int _\Omega w_i r_{m_i} d\Omega - \int _\Omega w_i  \frac{h_{ij}}{2} \frac{\partial r_{m_i}}{\partial x_j}d\Omega=0 </math>
938
|}
939
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
940
|}
941
942
where <math display="inline">w_i</math> are standard weighting functions.
943
944
Integrating by parts the second integral in ([[#eq-50|50]]) assuming that <math display="inline">w_i=0</math> on <math display="inline">\Gamma _v</math> and using the modified Neumann boundary conditions ([[#eq-49|49]]) we obtain
945
946
<span id="eq-51"></span>
947
{| class="formulaSCP" style="width: 100%; text-align: left;" 
948
|-
949
| 
950
{| style="text-align: left; margin:auto;width: 100%;" 
951
|-
952
| style="text-align: center;" | <math>\int _\Omega w_i r_{m_i} d\Omega + \sum \limits _{e=1}^{N_e} \int _{\Omega ^e} \frac{h_{ij}}{2} \frac{\partial w_i}{\partial x_j}r_{m_i}d\Omega=0 </math>
953
|}
954
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
955
|}
956
957
The characteristic distances can be defined in a number of ways. For instance choosing
958
959
<span id="eq-52"></span>
960
{| class="formulaSCP" style="width: 100%; text-align: left;" 
961
|-
962
| 
963
{| style="text-align: left; margin:auto;width: 100%;" 
964
|-
965
| style="text-align: center;" | <math>h_{ij}=h^e \frac{v_i}{\vert {v}\vert } </math>
966
|}
967
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
968
|}
969
970
where <math display="inline">h^e</math> is a characteristic distance for element <math display="inline">e</math>, reproduces the standard SUPG method <span id='citeF-11'></span><span id='citeF-39'></span>[[#cite-11|[11,39]]].
971
972
The effect of sharp internal gradients in the solution can be introduced by defining <math display="inline">h_{ij}</math> as
973
974
<span id="eq-53"></span>
975
{| class="formulaSCP" style="width: 100%; text-align: left;" 
976
|-
977
| 
978
{| style="text-align: left; margin:auto;width: 100%;" 
979
|-
980
| style="text-align: center;" | <math>h_{ij}=h^e_v\frac{v_i}{\vert {v}\vert } + h_g^e \frac{1}{\vert {v}\vert } {\partial v_{i} \over \partial x_j} </math>
981
|}
982
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
983
|}
984
985
The first term in the r.h.s. of Eq.([[#eq-53|53]]) introduces the SUPG stabilization while the second term introduces a shock-capturing type of stabilization. In Eq.([[#eq-53|53]]) <math display="inline">h^e_v</math> and <math display="inline">h_g^e</math> are appropriate characteristic distances for 1D elements <math display="inline">h_v^e = h_g^e = l^e</math>. The simplest choice for simplicial 2D and 3D elements is <math display="inline">h^e_v=h_g^e= 2 [\Omega ^e]^{1/n_s}</math>.
986
987
As for the treatment of the mass balance equation, the same procedure explained in Section 5.2 has been applied. Hence, the  stabilized mass balance expression is identical to Eq.([[#eq-19|19]]). Also, the corresponding  variational form of the mass balance equation coincides with Eq.([[#eq-35|35]]).
988
989
The final system of  discretized equations using  P1/P0+ elements is identical to Eq.([[#6 FEM DISCRETIZATION|6]]). For 3-noded triangles, all matrices and vectors coincide with the expressions of Box 1 except <math display="inline">{K}</math> and <math display="inline">{f}_v</math> which are now given by (for each element)
990
991
<span id="eq-54"></span>
992
{| class="formulaSCP" style="width: 100%; text-align: left;" 
993
|-
994
| 
995
{| style="text-align: left; margin:auto;width: 100%;" 
996
|-
997
| style="text-align: center;" | <math>{\bf K}^e ={\bf K}^e_0 + {\bf K}^e_v +  {\bf K}^e_s \quad ,\quad {\bf f}^e_v ={\bf f}^e_{v_0}  +  {\bf f}^e_s  </math>
998
|}
999
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
1000
|}
1001
1002
In Eq.([[#eq-54|54]]) <math display="inline">{\bf K}^e_0</math> is the viscous contribution  coinciding with the expression of <math display="inline">{\bf K}^e</math> of Box 1 and <math display="inline"> {\bf K}^e_v</math> and <math display="inline"> {\bf K}^e_s</math> are the convective and stabilization contributions for the element given by
1003
1004
<span id="eq-55.a"></span>
1005
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1006
|-
1007
| 
1008
{| style="text-align: left; margin:auto;width: 100%;" 
1009
|-
1010
| style="text-align: center;" | <math>{\mathbf{K}}_{v_{ij}}^e = \int _{\Omega ^e} \rho {N}_i ({v}^T {\boldsymbol \nabla } N_j){I}_{ns}  d\Omega \quad ,\quad  {K}^e_{s_{ij}}=  \int _{\Omega ^e} \frac{1}{2} \mathbf{N}_i \mathbf{G}^T \mathbf{H} (\mathbf{\nabla } {v}^T)^T \mathbf{N}_jd\Omega  </math>
1011
|}
1012
| style="width: 5px;text-align: right;white-space: nowrap;" | (55.a)
1013
|}
1014
1015
where <math display="inline">{\bf I}_{ns}</math> is the <math display="inline">n_s \times n_s</math> unit matrix and
1016
1017
<span id="eq-55.b"></span>
1018
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1019
|-
1020
| 
1021
{| style="text-align: left; margin:auto;width: 100%;" 
1022
|-
1023
| style="text-align: center;" | <math>\displaystyle{\bf G} = \left[\begin{matrix}{\boldsymbol \nabla } & {0} & {0}\\ {0}&{\boldsymbol \nabla } & {0}\\  {0} & {0}& {\boldsymbol \nabla } \end{matrix}  \right]~,~ \mathbf{H} = \left[\begin{matrix}{\bf h}_1 & {\bf 0}& {\bf 0}\\ {\bf 0}& {\bf h}_2 & {\bf 0}\\ {\bf 0} & {\bf 0}& {\bf h}_3 \end{matrix}  \right]~,~ \displaystyle {\bf h}_i = \left\{\begin{matrix}h_{i_1}\\h_{i_2}\\h_{i_3} \end{matrix}  \right\}~,~ {\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\} </math>
1024
|}
1025
| style="width: 5px;text-align: right;white-space: nowrap;" | (55.b)
1026
|}
1027
1028
In Eq.([[#eq-54|54]]) <math display="inline">{\bf f}^e_{v_0}</math> coincides with the expression of <math display="inline">{\bf f}^e_v</math> in Box 1 and <math display="inline"> {\bf f}^e_s</math> is the contribution from the stabilization terms given by
1029
1030
<span id="eq-56"></span>
1031
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1032
|-
1033
| 
1034
{| style="text-align: left; margin:auto;width: 100%;" 
1035
|-
1036
| style="text-align: center;" | <math>{\bf f}^e_{s_i} = \int _{\Omega ^e} \frac{1}{2} \mathbf{N}_i^T \mathbf{G}^T \mathbf{H} \mathbf{b} d\Omega   </math>
1037
|}
1038
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
1039
|}
1040
1041
The general form of <math display="inline">h_{ij}</math> of Eq.([[#eq-53|53]]) can be used in ([[#eq-56|56]]) to introduce SUPG and/or shock capturing effects in the stabilized formulation, as appropriate.
1042
1043
For the problems we are solving in this work the velocity field is free from internal layers. Therefore, good results are obtained with the standard Galerkin form of the momentum equations neglecting the stabilization term. The full stabilized formulation expressed by Eq.([[#eq-51|51]]), with <math display="inline">h_{ij}</math> given by Eq.([[#eq-53|53]]), should however be used for the general solution of Navier-Stokes problems involving high gradients of the velocity.
1044
1045
==9 PARTICULARIZATION FOR LAGRANGIAN FLOWS==
1046
1047
The formulation presented in Sections 2&#8211;6 for Stokes flows in an Eulerian framework can be easily particularized for analysis of a viscous fluid flow  using a Lagrangian description of the motion <span id='citeF-4'></span>[[#cite-4|[4]]].
1048
1049
The relevant change is the definition of the acceleration term in the momentum equations. In the ''updated Lagrangian'' formulation these are written as <span id='citeF-4'></span>[[#cite-4|[4]]]
1050
1051
<span id="eq-57"></span>
1052
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1053
|-
1054
| 
1055
{| style="text-align: left; margin:auto;width: 100%;" 
1056
|-
1057
| style="text-align: center;" | <math>r_{m_i} =\rho \frac{Dv_i}{Dt}- {\partial \sigma _{ij} \over \partial {}^{n+1} x_j}- {}^{n+1}b_i =0\quad i,j = 1,n_s </math>
1058
|}
1059
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
1060
|}
1061
1062
where <math display="inline">\frac{Dv_i}{Dt}</math> is the ''material derivative'' of the <math display="inline">i</math>th velocity component <math display="inline">v_i ({}^n {\bf x})</math> of a material point with coordinates <math display="inline">{}^n {\bf x}</math> at time <math display="inline">t=t_n</math>. Also, the Cauchy stresses <math display="inline">\sigma _{ij}</math>,  the body forces <math display="inline">b_i</math> and the coordinates <math display="inline">{}^{n+1}x_i</math> are referred to the ''updated configuration'' at time <math display="inline">t_{n+1}</math> <span id='citeF-4'></span>[[#cite-4|[4]]].
1063
1064
The material derivative of the velocity in the Lagrangian formulation is typically approximated as
1065
1066
<span id="eq-58"></span>
1067
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1068
|-
1069
| 
1070
{| style="text-align: left; margin:auto;width: 100%;" 
1071
|-
1072
| style="text-align: center;" | <math>\frac{Dv_i}{Dt} = \frac{{}^{n+1}v_i({}^{n+1}{\bf x})-{}^nv_i({}^n{\bf x})}{\Delta t} </math>
1073
|}
1074
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
1075
|}
1076
1077
where <math display="inline">{}^nv_i({}^n{\bf x})</math> denotes the value of the <math display="inline">i</math>th velocity component of the material point <math display="inline">{}^n {x}</math> at time <math display="inline">t_n</math>, etc.
1078
1079
'''Remark 10.''' The material derivative in an Eulerian description is expressed in a ''fixed control domain'' as
1080
1081
<span id="eq-59"></span>
1082
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1083
|-
1084
| 
1085
{| style="text-align: left; margin:auto;width: 100%;" 
1086
|-
1087
| style="text-align: center;" | <math>\frac{Dv_i}{Dt} =  {\partial v_i \over \partial t} +v_j  {\partial v_i \over \partial x_j}       </math>
1088
|}
1089
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
1090
|}
1091
1092
Eq.([[#eq-59|59]]) is the form used in the Eulerian formulation of the Navier-Stokes equations (see Eq.([[#eq-47|47]])). Neglecting the convective acceleration term in Eq.([[#eq-59|59]]) yields the expression used in Eq.([[#eq-1|1]]) for Stokes flows.
1093
1094
A difficulty in the Lagrangian formulation of flow problems is the need for tracking the motion of the material points in space and time. The authors have developed in recent years a particular class of Lagrangian method for analysis of fluid flow problems termed the Particle Finite Element Method (PFEM, www.cimne.com/pfem). The PFEM treats the nodes in a continuum as virtual particles that can freely move and even separate from the main domain representing, for instance, the effect of water drops in a splashing fluid or soil particles in an excavation problem. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved using the FEM. An advantage of the Lagrangian formulation in the PFEM is that the convective terms disappear from the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh node using a fast mesh regeneration procedure at each time step. The theory and applications of the PFEM are reported in <span id='citeF-1'></span><span id='citeF-8'></span>[[#cite-1|[1,8]]],<span id='citeF-18'></span>[[#cite-18|[18]]]&#8211;<span id='citeF-21'></span>[[#cite-21|[21]]],<span id='citeF-25'></span><span id='citeF-32'></span>[[#cite-25|[25,32]]],<span id='citeF-36'></span>[[#cite-36|[36]]]&#8211;<span id='citeF-38'></span>[[#cite-38|[38]]].
1095
1096
The particularization of the formulation of the P1/P0+ triangle presented in Sections 2&#8211;6 to the Lagrangian analysis of viscous flows is straightforward and simply implies computing the acceleration term with the expression given by Eq.([[#eq-58|58]]). The rest of the terms are identical to those given in those sections.
1097
1098
In this work we have applied the Lagrangian formulation of the P1/P0+ triangle to the study of multifluid problems with relatively small changes of the geometry of the fluid domain. The application of the P1/P0+ triangle and tetrahedra to a wider class of Lagrangian flow problems using the PFEM will be the subject of further work of the authors.
1099
1100
==10 EXAMPLES==
1101
1102
In this section we present the solutions obtained by the FIC method for four benchmark problems. The first two examples are steady state problems formulated in the Eulerian framework. The last two examples are transient problems formulated in the Lagrangian framework. In all these examples, the pressure solution has a strong discontinuity across the interface where the fluid viscosities are relevant and distinct. These pressure jumps occur either due to a relevant normal derivative of the velocity or in its absence due to a prescribed over-pressure force at the interface.
1103
1104
<div id='img-4'></div>
1105
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1106
|-
1107
|[[Image:Draft_Samper_536119838-fig_ex1_problem.png|294px|]]
1108
|[[Image:Draft_Samper_536119838-fig_ex4_problem.png|294px|]]
1109
|-style="text-align: center; font-size: 75%;"
1110
| (a) Eulerian extrusion problem. 
1111
| (b) Zhong's problem.
1112
|- style="text-align: center; font-size: 75%;"
1113
| colspan="2" | '''Figure 4:''' Description of the steady-state examples. (a) fixed-mesh (Eulerian) extrusion problem and (b) the Zhong's problem.
1114
|}
1115
1116
===10.1 Fixed-mesh (Eulerian) extrusion problem===
1117
1118
This example deals with a flow entering a rectangular 2D domain from the left boundary with a unit normal velocity and moving to the right where it finds an impermeable slip wall that deviates the flow upwards and downwards <span id='citeF-20'></span>[[#cite-20|[20]]]. The 2D domain occupies the region: <math display="inline">(x,y) \in [0, 1] \times [-0.5, 0.5]</math> (Figure 4a). The upper-half and lower-half of the domain are occupied by two immiscible fluids which have the same density (<math display="inline">\rho _{1} = \rho _{2} = \rho </math>), but with distinct viscosities (<math display="inline">\mu _{1} \neq \mu _{2}</math>). This problem has an analytical solution which can be expressed as follows.
1119
1120
<span id="eq-60"></span>
1121
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1122
|-
1123
| 
1124
{| style="text-align: left; margin:auto;width: 100%;" 
1125
|-
1126
| style="text-align: center;" | <math>{\bf v}(x,y) = \begin{bmatrix}1 - x \\  y \end{bmatrix} , \quad  p(x,y) = \left\{ \begin{array}{ll}\rho \left[x - \dfrac{1}{2}(x^{2} + y^{2}) -\dfrac{7}{24}\right]+ (\mu _{1} - \mu _{2}) & \forall  y > 0
1127
1128
\\  \rho \left[x - \dfrac{1}{2}(x^{2} + y^{2}) -\dfrac{7}{24}\right]- (\mu _{1} - \mu _{2}) & \forall  y < 0 \end{array} \right. </math>
1129
|}
1130
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
1131
|}
1132
1133
As the pressure solution is known up to a constant, in Eq. ([[#eq-60|60]]) the later is chosen such that the pressure has a zero mean value over the domain. At the interface (<math display="inline">y = 0</math>), due to a jump in the viscosities and a relevant directional derivative of the velocity along its normal, the pressure field has a jump across it.
1134
1135
'''Remark 11.''' This is a steady state example involving the solution of the Navier-Stokes equations solved in the Eulerian framework (Section 8). Generally, the presence of the convective term in the momentum equation triggers both global and local numerical instabilities which needs to be stabilized. However, these instabilities are manifested only when the solution of the continuous problem involves layers. In examples like the current one velocity field does not contain layers. Hence, in this problem, we can obtain a good numerical approximation to the velocity even without convective stabilization.
1136
1137
<div id='img-5a'></div>
1138
<div id='img-5'></div>
1139
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1140
|-
1141
|[[Image:Draft_Samper_536119838-fig_ex1_smesh_tri_grid_30_p_h.png|294px|]]
1142
|[[Image:Draft_Samper_536119838-fig_ex1_smesh_tri_grid_30_p_oL2.png|294px|]]
1143
|-style="text-align: center; font-size: 75%;"
1144
|(a) Obtained solution <math>p^h</math>; structured mesh.
1145
|(b) Best approximation <math>{\cal P}^h_{L^2}p</math>; structured mesh.
1146
|-
1147
|[[Image:Draft_Samper_536119838-fig_ex1_umesh_tri_grid_30_p_h.png|294px|]]
1148
|[[Image:Draft_Samper_536119838-fig_ex1_umesh_tri_grid_30_p_oL2.png|294px|]]
1149
|-style="text-align: center; font-size: 75%;"
1150
|(c) Obtained solution <math>p^h</math>; unstructured mesh.
1151
|(d) Best approximation <math>{\cal P}^h_{L^2}p</math>; unstructured mesh.
1152
|- style="text-align: center; font-size: 75%;"
1153
| colspan="2" | '''Figure 5:''' Fixed-mesh (Eulerian) extrusion problem. The obtained pressure solutions compared with the corresponding best approximations. The material properties are chosen as <math>(\rho _{1}, \mu _{1}) = (5,5)</math> and <math>(\rho _{2}, \mu _{2}) = (5,1)</math>.
1154
|}
1155
1156
Figure 5a illustrates the pressure solution <math display="inline">p^{h}</math> obtained by the P1-P0+ triangle. The domain is discretized using a symmetric structured mesh of <math display="inline">2 \times 30 \times 30</math> elements. The material properties are chosen to be: <math display="inline">\rho = 5</math>, <math display="inline">\mu _{1} = 5</math> and <math display="inline">\mu _{2} = 1</math>. The velocity boundary conditions are chosen to be of Dirichlet (essential) type, taking values given by the expression in Eq. ([[#eq-60|60]]). The indeterminacy of the pressure solution is removed by imposing the zero-mean pressure condition <span id='citeF-11'></span>[[#cite-11|[11]]]. The nonlinear convective terms in the Navier-Stokes equations are linearized using the Newton&#8211;Raphson method (Appendix A) and they converge in just two iterations for a tolerance of <math display="inline">10^{-6}</math>.
1157
1158
Figure 5b illustrates the best approximation to the exact pressure solution <math display="inline">p</math> given in Eq. ([[#eq-60|60]]) from the discrete pressure space <math display="inline">Q^{h}</math> of piecewise constant pressure values over the elements of the considered mesh. The best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math> is chosen from <math display="inline">Q^{h}</math> using the <math display="inline">\hbox{L}^{2}</math> norm (<math display="inline">\| \star \| := \sqrt{\int _{\Omega } (\star )^{2} \mbox{d}\Omega }</math>) as the metric. Thus, <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math> is the <math display="inline">\hbox{L}^{2}</math> projection of <math display="inline">p</math> onto <math display="inline">Q^{h}</math>. It can be obtained as follows: Find <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p \in Q^{h}</math>, such that
1159
1160
<span id="eq-61"></span>
1161
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1162
|-
1163
| 
1164
{| style="text-align: left; margin:auto;width: 100%;" 
1165
|-
1166
| style="text-align: center;" | <math>\int _{\Omega }{q^{h}(p-\mathcal{P}^{h}_{\hbox{L}^{2}}p) \mbox{d}\Omega } = 0 \forall  q^{h} \in Q^{h} \Rightarrow  {\|{p-\mathcal{P}^{h}_{\hbox{L}^{2}}p}\| } \leq {\|{p-p^{h}}\| } \forall  p^{h} \in Q^{h} </math>
1167
|}
1168
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
1169
|}
1170
1171
Figure 5c illustrates the constant pressure solution <math display="inline">p^{h}</math> (termed <math display="inline">\bar p^e</math> in Section 6) obtained by the new  P1-P0+ triangle and using an unstructured mesh. The unstructured mesh is generated by small random perturbations of the interior nodes (excluding the nodes on the interface) of the structured mesh followed by a Delaunay tessellation. Figure [[#img-5a|5d]] illustrates the best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math> obtained using the same mesh.
1172
1173
<div id='img-6'></div>
1174
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1175
|-
1176
|[[Image:Draft_Samper_536119838-fig_ex1_smesh_tri_grids_L2_vol_err.png|294px|]]
1177
|[[Image:Draft_Samper_536119838-fig_ex4_smesh_tri_grids_L2_vol_err.png|294px|]]
1178
|-style="text-align: center; font-size: 75%;"
1179
|(a) Fixed-mesh (Eulerian) extrusion problem.
1180
|(b) Zhong's problem.
1181
|- style="text-align: center; font-size: 75%;"
1182
| colspan="2" | '''Figure 6:''' Error convergence in the <math>\hbox{L}^{2}</math> norm on mesh refinement for <math>p^{h}</math> and <math>(\boldsymbol \nabla  \cdot {v}^{h})</math>.
1183
|}
1184
1185
Figure 6a shows the convergence of the error <math display="inline">(p - p^{h})</math> in the <math display="inline">\hbox{L}^{2}</math> norm with respect to mesh refinement. This error is compared with the best approximation error <math display="inline">(p - \mathcal{P}^{h}_{\hbox{L}^{2}}p)</math>. To measure the convergence rate the domain is discretized by a sequence of symmetrical structured meshes of <math display="inline">2 \times n \times n</math> elements, using <math display="inline">n \in \{ 20,24,28,32,36,40,44,48,52,56\} </math>. The relative error in the pressure solutions are measured as follows.
1186
1187
<span id="eq-62"></span>
1188
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1189
|-
1190
| 
1191
{| style="text-align: left; margin:auto;width: 100%;" 
1192
|-
1193
| style="text-align: center;" | <math>e^{h}_{p\hbox{L}^{2}} = \frac{\| p - p^{h}\| }{\| p\| } = \frac{\sqrt{\int _{\Omega } (p-p^{h})^{2} \mbox{d}\Omega }}{\sqrt{\int _{\Omega } p^{2} \mbox{d}\Omega }} </math>
1194
|}
1195
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
1196
|}
1197
1198
As expected for a piecewise constant pressure approximation, the convergence rate is found to be one for both solutions. Further, the error line of the numerical solution <math display="inline">p^{h}</math> is found to be close to the error line of the best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math>, indicating the good accuracy of the obtained solution. Additionally, in Figure 6a the convergence in the <math display="inline">\hbox{L}^{2}</math> norm of the divergence of the velocity is also shown. As the exact solution is solenoidal, this error is measured as follows.
1199
1200
<span id="eq-63"></span>
1201
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1202
|-
1203
| 
1204
{| style="text-align: left; margin:auto;width: 100%;" 
1205
|-
1206
| style="text-align: center;" | <math>\| \boldsymbol \nabla  \cdot {\bf v}^{h} \| = \sqrt{\int _{\Omega } (\boldsymbol \nabla  \cdot {\bf v}^{h})^{2} \mbox{d}\Omega } </math>
1207
|}
1208
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
1209
|}
1210
1211
where <math display="inline">{\bf v}^{h}</math> is the velocity solution obtained with the P1/P0+ triangle. The convergence rate for <math display="inline">\| \boldsymbol \nabla  \cdot {\bf v}^{h} \| </math> is found to be close to <math display="inline">1.5</math>, which is more than the expected first-order rate. This super-convergence for <math display="inline">\| \boldsymbol \nabla  \cdot {\bf v}^{h} \| </math> could be due to the smooth (linear, cf. Eq. ([[#eq-60|60]])) solution profile for the velocity.
1212
1213
===10.2 Zhong's problem: buoyant Stokes flow with a columnar viscosity structure===
1214
1215
This example deals with the buoyancy-driven Stokes flow of a fluid with a columnar viscosity structure and was first presented in <span id='citeF-40'></span>[[#cite-40|[40]]]. It has been previously used as a benchmark in <span id='citeF-13'></span><span id='citeF-24'></span>[[#cite-13|[13,24]]]. It is an idealization of the mantle convection and the plate dynamics that occur in the subduction zones of the Earth. The temperature of the subducted plates are much colder than the ambient mantle which in turn leads to sharp lateral variations in the viscosity. The 2D domain (cf. Figure [[#img-4|4b]]) is a unit square: <math display="inline">(x,y) \in [0, 1] \times [0, 1]</math>. The left-half and right-half of the domain have viscosities given by <math display="inline">\mu _{1} = 1</math> and <math display="inline">\mu _{2} = 10^6</math>, respectively. The density and the acceleration due to gravity are specified as <math display="inline">\rho (x,y) = - \cos (\pi x) \sin (\pi y)</math> and <math display="inline">{\bf g} = (0,-1)</math>, respectively. Free-slip boundary conditions are imposed everywhere on the domain boundary.
1216
1217
The evaluation <span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]] of the analytical solution is described in <span id='citeF-40'></span>[[#cite-40|[40]]] using the method of separation of variables (used earlier in <span id='citeF-16'></span>[[#cite-16|[16]]]) and the propagator matrix techniques <span id='citeF-15'></span>[[#cite-15|[15]]]. The constant to fix the pressure is chosen such that the latter has a zero mean value over the domain. Due to a jump in the viscosities at the viscosity boundary (<math display="inline">x = 0.5</math>) and a non-zero directional derivative of the velocity along its normal, the pressure field acquires a jump across it. The expression for the pressure jump at the viscosity boundary can be written as follows.
1218
1219
<span id="eq-64"></span>
1220
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1221
|-
1222
| 
1223
{| style="text-align: left; margin:auto;width: 100%;" 
1224
|-
1225
| style="text-align: center;" | <math>[[p (x=0.5,y)  ]] = \frac{2 |{g}| (\mu _{2}-\mu _{1})^2 \cosh ^2(\pi /2) \cos (\pi y)}{(\mu _{1}+\mu _{2})^2 \sinh ^2(\pi ) - \pi ^2 (\mu _{2}-\mu _{1})^2} </math>
1226
|}
1227
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
1228
|}
1229
1230
<div id='img-7'></div>
1231
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1232
|-
1233
|[[Image:Draft_Samper_536119838-fig_ex4_smesh_tri_grid_30_p_h.png|294px|]]
1234
|[[Image:Draft_Samper_536119838-fig_ex4_smesh_tri_grid_30_p_oL2.png|294px|]]
1235
|-style="text-align: center; font-size: 75%;"
1236
|(a) Obtained solution <math>p^h</math>; structured mesh.
1237
|(b) Best approximation <math>{\cal P}^h_{L^2}p</math>; structured mesh.
1238
|-
1239
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_30_p_h.png|294px|]]
1240
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_30_p_oL2.png|294px|]]
1241
|-style="text-align: center; font-size: 75%;"
1242
|(c) Obtained solution <math>p^h</math>; unstructured mesh.
1243
|(d) Best approximation <math>{\cal P}^h_{L^2}p</math>; unstructured mesh.
1244
|- style="text-align: center; font-size: 75%;"
1245
| colspan="2" | '''Figure 7:''' Zhong's problem: buoyancy-driven Stokes flow with a columnar viscosity structure in a square domain with slip walls. The density, the acceleration due to gravity and the viscosities are specified as <math>\rho = -\cos (\pi x) \sin (\pi y)</math>, <math>{g} = (0,-1)</math>, <math>\mu _{1} = 1</math> and <math>\mu _{2} = 10^6</math>, respectively.
1246
|}
1247
1248
Figures 7a and 7c show the pressure solution <math display="inline">p^h</math> obtained by the P1-P0+ triangle using a structured mesh (<math display="inline">2 \times 30 \times 30</math> elements) and an unstructured mesh (<math display="inline">2 \times 30 \times 30</math> elements), respectively. The corresponding best approximations in the <math display="inline">\hbox{L}^{2}</math> norm <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math> are shown in figures  fig_ex4_smesh_tri_grid_30_p_oL2 and [[#img-7|7]], respectively.
1249
1250
Figure [[#img-6|6b]] shows the convergence of the error <math display="inline">(p - p^{h})</math> in the <math display="inline">\hbox{L}^{2}</math> norm with respect to mesh refinement. This error is compared with the best approximation error <math display="inline">(p - \mathcal{P}^{h}_{\hbox{L}^{2}}p)</math>. As expected for a piecewise constant pressure approximation, the convergence rate is found to be one for both solutions. Further, the error line for the numerical solution <math display="inline">p^{h}</math> is found to be close to the error line of the best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p</math>, indicating the good accuracy of the obtained solution. Additionally, in Figure [[#img-6|6b]] the convergence in the <math display="inline">\hbox{L}^{2}</math> norm of the divergence of the velocity is also shown.  As expected for a piecewise linear velocity approximation, a first-order convergence rate is found for <math display="inline">\| \boldsymbol \nabla  \cdot {\bf v}^{h} \| </math>.
1251
1252
<span id="fn-1"></span>
1253
<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) The algebraic work involved to arrive at a closed form analytical expression of the solution is very tedious which includes the symbolic inversion of a <math>4 \times 4</math> matrix. Thus, the analytical solution is often computed numerically for a specified problem data.</span>
1254
1255
===10.3 Moving-mesh (Lagrangian) extrusion problem===
1256
1257
This is a transient example where the initial configuration is a rectangular 2D domain occupied by two immiscible fluids contained (due to the action of gravity) in a region bounded with sufficiently high walls (Figure 8a). Both the walls and the floor are considered to be slippery. The floor and the left wall are fixed while the right wall is pushed towards left at a constant prescribed velocity <math display="inline">V_{p}</math>. This is a free surface problem with an interface between the two immiscible fluids and the shape of the domain varies due to the action of the right wall. The numerical solution is sought using a Lagrangian formulation of the problem (Section 9).
1258
1259
<div id='img-8a'></div>
1260
<div id='img-8'></div>
1261
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1262
|-
1263
|[[Image:Draft_Samper_536119838-fig_ex2_problem.png|294px|]]
1264
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_initial.png|294px|]]
1265
|-style="text-align: center; font-size: 75%;"
1266
| (a) Lagrangian) extrusion problem.
1267
| (b) Lagrangian example with an unstable interface.
1268
|- style="text-align: center; font-size: 75%;"
1269
| colspan="2" | '''Figure 8:''' Description of the transient examples. (a) moving-mesh (Lagrangian) extrusion problem and (b) Lagrangian example with fixed-walls and an unstable interface: the interface initially has a serrated profile.
1270
|}
1271
1272
<div id='img-9a'></div>
1273
<div id='img-9'></div>
1274
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1275
|-
1276
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_initial.png|294px|]]
1277
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_final.png|294px|]]
1278
|-style="text-align: center; font-size: 75%;"
1279
| (a)  initial state (t = 0); structured mesh.
1280
|(b) final state (t = 0); structured mesh.
1281
|-
1282
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_initial.png|294px|]]
1283
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_final.png|294px|]]
1284
|-style="text-align: center; font-size: 75%;"
1285
|(c) initial state (t = 0); unstructured mesh.
1286
|(d) final state (t = 0); unstructured mesh.
1287
|- style="text-align: center; font-size: 75%;"
1288
| colspan="2" | '''Figure 9:''' Moving-mesh (Lagrangian) extrusion problem: the initial (<math>t = 0</math>) and final (<math>t = 2</math>) configurations are independent to the choice of the mesh-type.
1289
|}
1290
1291
The fluid viscosities are chosen as <math display="inline">\mu _{1} = 1</math> and <math display="inline">\mu _{2} = 10</math>. The fluid densities are chosen as <math display="inline">\rho _{1} = 1</math> and <math display="inline">\rho _{2} = 5</math>. The acceleration due to gravity is taken as <math display="inline">[0,-10] \hbox{m/s}^{2}</math> and the prescribed velocity <math display="inline">V_{p} = 0.1 \hbox{m/s}</math>. The 2D domain has an initial length <math display="inline">l_{1} = 0.8 \hbox{m}</math> and height <math display="inline">l_{2} = 0.4 \hbox{m}</math>. The upper-half of the domain is occupied by the fluid with properties <math display="inline">\rho _{1}, \mu _{1}</math>. For these data, the dynamics is dominated by the viscous term and hence in the exact solution the inertial effects can be neglected.
1292
1293
As a shorthand notation, we denote by <math display="inline">{}^{0}{\bf x}</math> and <math display="inline">{}^{t}{\bf x}</math> as the positions of an arbitrary fluid particle at times <math display="inline">t = 0</math> and <math display="inline">t</math>, respectively. Similar notation is used for the rest of the dependent variables in the Lagrangian description, e.g. <math display="inline">{}^{0}{\bf v}</math> and <math display="inline">^{t}{\bf v}</math>. Using this notation and the expressions <math display="inline">L_{1} := (l_{1} - V_{p}t)</math> and <math display="inline">L_{2} := (l_{1}l_{2}/L_{1})</math>,  the ''exact solution'' to this problem can be expressed as follows.
1294
1295
<span id="eq-65"></span>
1296
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1297
|-
1298
| 
1299
{| style="text-align: left; margin:auto;width: 100%;" 
1300
|-
1301
| style="text-align: center;" | <math>{}^{t}{\bf x} = \begin{bmatrix}(L_{1}/l_{1}) & 0 \\  0 & (L_{2}/l_{2}) \end{bmatrix} {}^{0}{\bf v},\quad  {}^{t}{\bf v} := \frac{\mbox{d}}{\mbox{d}t}({\bf x}^{t}) = \frac{v_{p}}{L_{1}} \begin{bmatrix}-1 & 0 \\  0 & 1 \end{bmatrix} {}^{t}{\bf v} </math>
1302
|}
1303
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
1304
|}
1305
1306
<span id="eq-66"></span>
1307
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1308
|-
1309
| 
1310
{| style="text-align: left; margin:auto;width: 100%;" 
1311
|-
1312
| style="text-align: center;" | <math>{}^{t}{\bf a} := \frac{\mbox{d}}{\mbox{d}t}({}^{t}{\bf v}) = \left(\frac{V_{p}}{L_{1}}\right)^{2} \begin{bmatrix}0 & 0 \\  0 & 2 \end{bmatrix} {}^{t}{\bf x},\quad  {}^{t}\boldsymbol \nabla ({}^{t}{\bf v}) := \frac{\mbox{d}}{\mbox{d}{}^{t}{\bf x}}({}^{t}{\bf v}) = \frac{V_{p}}{L_{1}} \begin{bmatrix}-1 & 0 \\  0  & 1 \end{bmatrix} </math>
1313
|}
1314
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
1315
|}
1316
1317
<span id="eq-67"></span>
1318
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1319
|-
1320
| 
1321
{| style="text-align: left; margin:auto;width: 100%;" 
1322
|-
1323
| style="text-align: center;" | <math>p ({\bf x},t) = \left\{ \begin{array}{ll}\rho _{1} g (L_{2}-{}^{t}Y) & \hbox{ for} L_{2} \geq ^{t}Y > \dfrac{L_{2}}{2}
1324
1325
\\  \rho _{1} g\dfrac{L_{2}}{2} + 2(\mu _{2}-\mu _{1})\dfrac{V_{p}}{L_{1}} + \rho _{2}g\left(\dfrac{L_{2}}{2} - ^{t}Y\right)& \hbox{ for} \dfrac{L_{2}}{2} > ^{t}Y \geq 0 \end{array} \right. </math>
1326
|}
1327
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
1328
|}
1329
1330
where <math display="inline">^{t}{\bf a}</math> denotes the acceleration of the fluid particle and <math display="inline">^{t}Y</math> is the y-coordinate of the trajectory <math display="inline">^{t}{x}</math>. The constant in the pressure solution given by Eq. ([[#eq-67|67]]) is chosen such that it takes a value zero at the free surface. In the numerical solution, no boundary conditions are imposed for the velocity  at the free surface. On the rest of the boundary, Dirichlet (essential) conditions are imposed for the velocity. As the velocity boundary conditions are not exclusively of the Dirichlet type, no further conditions need to be imposed for the pressure. For the sake of simplicity the discretized (Lagrangian) equations are linearized using the Picard (Fixed-point) method. The configuration <math display="inline">^{t}{x}</math> is updated at every iteration. The time integration is done using the Newmark method with the choice <math display="inline">\theta = 1/2</math> and <math display="inline">\beta = 1/4</math> (zero dissipation; second-order accuracy). The time increment is chosen as <math display="inline">\Delta t = 0.1</math>s (Appendix B).
1331
1332
Figure 9b illustrates the final configuration (at time <math display="inline">t = 2 \hbox{s}</math>) obtained with the P1-P0+ triangle  using a symmetric structured mesh of <math display="inline">2 \times 24 \times 12</math> elements. The initial configuration is shown in Figure 9a. Figures 9c and [[#img-9a|9d]] show the initial and final configurations when using an unstructured mesh. For the initial domain area to be conserved the interface and the free surface should raise to a height of <math display="inline">0.2667 \hbox{m}</math> and <math display="inline">0.5333 \hbox{m}</math>, respectively. This is indeed observed in Figures 9b and [[#img-9a|9d]] showing the good area conservation properties for the P1/P0+ triangle.
1333
1334
<div id='img-10a'></div>
1335
<div id='img-10'></div>
1336
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1337
|-
1338
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_160_20.png|294px|]]
1339
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_160_20_BA.png|294px|]]
1340
|-style="text-align: center; font-size: 75%;"
1341
|(a) Obtained solution <math>p^{2,h}</math>; structured mesh.
1342
|(b) Best approximation <math>{\cal P}^h_{L^2}p^2</math>; structured mesh.
1343
|-
1344
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_160_20.png|294px|]]
1345
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_160_20_BA.png|294px|]]
1346
|-style="text-align: center; font-size: 75%;"
1347
|(c) Obtained solution <math>p^{2,h}</math>; unstructured mesh.
1348
|(d) Best approximation <math>{\cal P}^h_{L^2}p^2</math>; unstructured mesh.
1349
|- style="text-align: center; font-size: 75%;"
1350
| colspan="2" | '''Figure 10:''' Moving-mesh (Lagrangian) extrusion problem. The obtained pressure solutions at time <math>t = 2</math> are compared with the corresponding best approximations. The material properties are chosen as <math>(\rho _{1}, \mu _{1}) = (1,1)</math> and <math>(\rho _{2}, \mu _{2}) = (5,10)</math>.
1351
|}
1352
1353
Figure 10a shows the pressure solution <math display="inline">p^{2,h}</math> obtained using the structured mesh. As the two fluids have distinct densities, it leads to a jump in both the pressure and its gradient at the interface. Figure 10b shows the <math display="inline">\hbox{L}^{2}</math> projection of the exact pressure <math display="inline">p^{2}</math> given in Eq. ([[#eq-67|67]]) onto the discrete pressure space <math display="inline">Q^{h}</math> and is denoted as <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math>. For the sake of comparison, the constants in the discrete pressure solutions <math display="inline">p^{2,h}</math> and <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math> are chosen ''a posteriori'' such that the minimum value on the free surface is zero. Likewise, Figures 10c and [[#img-10a|10d]] show the pressure solutions <math display="inline">p^{2,h}</math> and <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math> obtained on  the unstructured mesh.
1354
1355
<div id='img-11a'></div>
1356
<div id='img-11'></div>
1357
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1358
|-
1359
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_0_0.png|294px|]]
1360
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_0_0_BA.png|294px|]]
1361
|-style="text-align: center; font-size: 75%;"
1362
|(a) Obtained solution <math>p^{2,h}</math>; structured mesh.
1363
|(b) Best approximation <math>{\cal P}^h_{L^2}p^2</math>; structured mesh.
1364
|-
1365
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_0_0.png|294px|]]
1366
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_0_0_BA.png|294px|]]
1367
|-style="text-align: center; font-size: 75%;"
1368
|(c) Obtained solution <math>p^{2,h}</math>; unstructured mesh.
1369
|(d) Best approximation <math>{\cal P}^h_{L^2}p^2</math>; unstructured mesh.
1370
|- style="text-align: center; font-size: 75%;"
1371
| colspan="2" | '''Figure 11:''' Moving-mesh (Lagrangian) extrusion problem: <math>p^{2,h}</math> and <math>\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math> projected onto the <math>x-z</math> plane. The material properties are chosen as <math>(\rho _{1}, \mu _{1}) = (1,1)</math> and <math>(\rho _{2}, \mu _{2}) = (5,10)</math>.
1372
|}
1373
1374
Figure [[#img-11|11]] presents the same results as shown in Figure [[#img-10|10]] but with the solutions projected onto the <math display="inline">x-z</math> plane. In this view, the inter-element jumps in the obtained solution <math display="inline">p^{2,h}</math> can be better compared with those in the best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math>. We see that on both meshes the obtained solution <math display="inline">p^{2,h}</math> is in good agreement with the best approximation <math display="inline">\mathcal{P}^{h}_{\hbox{L}^{2}}p^{2}</math>.
1375
1376
<div id='img-12a'></div>
1377
<div id='img-12'></div>
1378
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1379
|-
1380
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_160_20_enri.png|294px|]]
1381
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_160_20_enri.png|294px|]]
1382
|-style="text-align: center; font-size: 75%;"
1383
|(a) Effective pressure <math>p^{2,h}</math>; structured mesh.
1384
|(b) Effective pressure <math>p^{2,h}</math>; unstructured mesh.
1385
|- style="text-align: center; font-size: 75%;"
1386
| colspan="2" | '''Figure 12:''' Moving-mesh (Lagrangian) extrusion problem: The effective pressure solutions at time <math>t = 2</math> which includes a prescribed linear field defined over each element. The material properties are chosen as <math>(\rho _{1}, \mu _{1}) = (1,1)</math> and <math>(\rho _{2}, \mu _{2}) = (5,10)</math>.
1387
|}
1388
1389
Figures 12a and [[#img-12a|12b]] illustrate the ''effective'' pressure solution obtained with the P1/P0+ triangle using the structured and the unstructured meshes, respectively. Recall that the effective pressure solution includes a prescribed linear field defined over each element given by Eq.([[#eq-37|37]]). A nearly identical figure is obtained by plotting the exact pressure solution given by Eq.([[#eq-67|67]]), thus indicating the good accuracy of the method.
1390
1391
Finally, in Figure [[#img-13|13]] we present the time evolution of the error in the area for all the cases discussed in this section. The relative error in the area is measured in three different ways
1392
1393
<span id="eq-68"></span>
1394
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1395
|-
1396
| 
1397
{| style="text-align: left; margin:auto;width: 100%;" 
1398
|-
1399
| style="text-align: center;" | <math>e^{h}_{\Omega } := \dfrac{\displaystyle \sum _{e \in \Omega }({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e \in \Omega } {}^{0}\Omega _{e}} ,\quad  e^{h}_{\Omega _{1}} := \dfrac{\displaystyle \sum _{e \in \Omega _{1}}({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e \in \Omega _{1}} {}^{0}\Omega _{e}} ,\quad  e^{h}_{\Omega _{2}} := \dfrac{\displaystyle \sum _{e \in \Omega _{2}}({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e \in \Omega _{2}} {}^{0}\Omega _{e}} </math>
1400
|}
1401
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
1402
|}
1403
1404
where <math display="inline">{}^{0}\Omega _{e}</math> and <math display="inline">{}^t \Omega _{e}</math> represent the initial and current (at time <math display="inline">t</math>) areas of an element with index <math display="inline">e</math>. <math display="inline">\Omega _{1}</math> and <math display="inline">\Omega _{2}</math> denote the regions occupied by the two immiscible fluids, respectively and <math display="inline">\Omega :=  \Omega _{1} \cup \Omega _{2}</math> is the domain. The indices <math display="inline">1</math> and <math display="inline">2</math> refer to the upper-half and lower-half regions of the domain. <math display="inline">e^{h}_{\Omega }</math> is the relative error in the total area of both fluids. <math display="inline">e^{h}_{\Omega _{1}}</math> and <math display="inline">e^{h}_{\Omega _{2}}</math> are the relative errors in the total area of the fluids occupying the regions <math display="inline">\Omega _{1}</math> and <math display="inline">\Omega _{2}</math>, respectively. For the simulated time period (<math display="inline">2</math>s), the maximum relative error is below <math display="inline">3.25 \times 10^{-4}</math> for <math display="inline">e^{h}_{\Omega }</math>, <math display="inline">e^{h}_{\Omega _{1}}</math> and <math display="inline">e^{h}_{\Omega _{2}}</math>.
1405
1406
<div id='img-13a'></div>
1407
<div id='img-13'></div>
1408
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1409
|-
1410
|[[Image:Draft_Samper_536119838-fig_ex2_smesh_tri_grid_24_12_rho_diff_vol_err.png|294px|]]
1411
|[[Image:Draft_Samper_536119838-fig_ex2_umesh_tri_grid_24_12_rho_diff_vol_err.png|294px|]]
1412
|-style="text-align: center; font-size: 75%;"
1413
|(a) Area errors using a structured mesh.
1414
|(b) Area errors using an unstructured mesh.
1415
|- style="text-align: center; font-size: 75%;"
1416
| colspan="2" | '''Figure 13:''' Moving-mesh (Lagrangian) extrusion problem: time evolution of the area errors. The material properties are chosen as <math>(\rho _{1}, \mu _{1}) = (1,1)</math> and <math>(\rho _{2}, \mu _{2}) = (5,10)</math>.
1417
|}
1418
1419
===10.4 Lagrangian example with fixed-walls and an unstable interface===
1420
1421
Finally, we present another transient example where the two-fluid interface initially has a serrated profile (thus unstable) and becomes straight in due course. This configuration (Figure [[#img-8a|8b]]) is obtained by meshing the domain with an unstructured mesh and identifying the elements that occupy the upper-half of the domain to represent a fluid with lower density and the remaining elements to represent the other fluid with a higher density.
1422
1423
To observe the interface evolution in a reasonable time period, fluid-pairs with low viscosities are chosen. Further, to avoid sloshing-like effects that arise due to the low viscosities, all the walls of the domain are fixed and are considered slippery. Thus, there are only buoyancy driven low-speed motions in the fluid due to its initial unstable configuration.
1424
1425
Due to the small viscosities and low-speed motions, the pressure has a negligible jump across the interface. In order to have a visible pressure jump, we prescribe a fictitious over-pressure force at the interface given by the expression <math display="inline">{t}^{\hbox{I}} = -\gamma {n}</math>. Here <math display="inline">{n}</math> represents the outward normal to the interface with respect to the heavier fluid and <math display="inline">\gamma </math> represents the prescribed jump in the pressure value across the interface.
1426
1427
This problem is solved considering two fluid-pairs. For the first fluid-pair, the material properties of the fluids on top and bottom are taken as <math display="inline">(\rho _{1},\mu _{1}) = (1, 10^{-6})</math> and <math display="inline">(\rho _{2},\mu _{2}) = (5, 10^{-5})</math>, respectively. At the interface, the over-pressure force is taken as <math display="inline">\gamma = 5</math>. The second fluid-pair consists of sunflower oil on the top and water at the bottom. The material properties for sunflower oil and water at <math display="inline">25^{\circ } \hbox{C}</math> are (in International units) <math display="inline">(\rho _{1},\mu _{1}) = (919, 4.9 \times 10^{-2})</math> and <math display="inline">(\rho _{2},\mu _{2}) = (997, 8.9 \times 10^{-4})</math>, respectively. At the oil&#8211;water interface, the over-pressure force is taken as <math display="inline">\gamma = 1000</math>. The initial configuration is the same for both fluid-pairs, cf. figure [[#img-8a|8b]]. The time integration is done by the Newmark method with the choice <math display="inline">\theta = 1/2</math> and <math display="inline">\beta = 1/4</math> (zero-dissipation; second-order accuracy). The time increment is chosen as <math display="inline">\Delta t = 0.01</math>s. As no re-meshing is done, the simulation is stopped after <math display="inline">100</math> time steps, i.e. a simulation time of <math display="inline">1</math>s, to avoid numerical difficulties associated with severe mesh deformations.
1428
1429
<div id='img-14a'></div>
1430
<div id='img-14'></div>
1431
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1432
|-
1433
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_tl.png|294px|]]
1434
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_tl.png|294px|]]
1435
|-style="text-align: center; font-size: 75%;"
1436
|(a) First fluid-pair interface length evolution.
1437
|(b) Oil-Water interface length evolution.
1438
|-
1439
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_final.png|294px|]]
1440
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_final.png|294px|]]
1441
|-style="text-align: center; font-size: 75%;"
1442
|(c) First fluid-pair domain at ''t'' = 0.11s.
1443
|(d) Oil-Water domain at ''t'' = 0.55s.
1444
|- style="text-align: center; font-size: 75%;"
1445
| colspan="2" | '''Figure 14:''' Lagrangian example with fixed-walls and an unstable interface. The time evolution of the interface length and the domain when the former first attains a minimum are shown. The left column corresponds to a fluid-pair: <math>(\rho _{1},\mu _{1}) = (1, 10^{-6})</math>, <math>(\rho _{2},\mu _{2}) = (5, 10^{-5})</math> and <math>\gamma = 5</math>. The right column corresponds to oil&#8211;water interface: <math>(\rho _{1},\mu _{1}) = (919, 4.9 \times 10^{-2})</math> and <math>(\rho _{2},\mu _{2}) = (997, 8.9 \times 10^{-4})</math> and <math>\gamma = 1000</math>.
1446
|}
1447
1448
<div id='img-15a'></div>
1449
<div id='img-15'></div>
1450
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1451
|-
1452
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_160_20.png|294px|]]
1453
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_0_0.png|294px|]]
1454
|-style="text-align: center; font-size: 75%;"
1455
|(a) Obtained solution <math>P^{0.11,h}</math>; 3D view.
1456
|(b) <math>P^{0.11,h}</math> projected onto the <math>x - z</math> plane.
1457
|-
1458
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_160_20.png|294px|]]
1459
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_0_0.png|294px|]]
1460
|-style="text-align: center; font-size: 75%;"
1461
|(c) Obtained solution <math>P^{0.55,h}</math>; 3D view.
1462
|(d) <math>P^{0.55,h}</math> projected onto the <math>x - z</math> plane.
1463
|- style="text-align: center; font-size: 75%;"
1464
| colspan="2" | '''Figure 15:''' Lagrangian example with fixed-walls and an unstable interface. The obtained pressure solutions when the interface lengths first attain a minimum are shown. The top row corresponds to a fluid-pair: <math>(\rho _{1},\mu _{1}) = (1, 10^{-6})</math>, <math>(\rho _{2},\mu _{2}) = (5, 10^{-5})</math> and <math>\gamma = 5</math>. The bottom row corresponds to oil&#8211;water interface: <math>(\rho _{1},\mu _{1}) = (919, 4.9 \times 10^{-2})</math> and <math>(\rho _{2},\mu _{2}) = (997, 8.9 \times 10^{-4})</math> and <math>\gamma = 1000</math>.
1465
|}
1466
1467
Due to the small viscosities, the profile of the interface will have an oscillatory behavior (possibly damped due to momentum distribution in all directions) and a steady state can be expected in due course. We use the total length of the interface as a metric to measure the “distance” of the perturbed state of the interface from the stable horizontal profile. For an oscillating interface profile, the interface length will have an oscillatory behavior and will take a minimum value equal to <math display="inline">0.8</math> (horizontal length of the domain) whenever the interface profile is horizontal.
1468
1469
Within the simulated time period, the  P1-P0+ triangle is able to reproduce this oscillatory behavior of the interface for the first fluid-pair, cf. Figure 14a. The interface length first attains a minimum value of <math display="inline">0.80014</math> in eleven time steps (at <math display="inline">t = 0.11</math>s). The interface at this instant attains a horizontal profile as shown in Figure 14c. However, near the end of the simulation we observe an unphysical increase in the interface length. This could be due to the numerical difficulties associated to mesh-distortion (recall that no re-meshing is done here).
1470
1471
As the viscosities in the second fluid-pair (oil&#8211;water) are three orders of magnitude greater than those for the first fluid-pair, it is reasonable to expect a greater time scale for the former. Figure 14b shows that the oil&#8211;water interface length gradually reduces (i.e. the interface flattens) and takes a minimum value of <math display="inline">0.80842</math> at <math display="inline">t = 0.55</math>s. The interface at this instant attains a nearly horizontal profile as shown in Figure [[#img-14a|14d]]. Past <math display="inline">t = 0.55</math>s, the interface length gradually returns back to its initial value. However, the trend near the end of the simulation shows an unphysical increase in the interface length, which again could be due to the numerical difficulties associated to mesh-distortion.
1472
1473
Figure [[#img-15|15]] shows the pressure solutions obtained for the considered fluid-pairs when their interfaces first attain a minimum length. We observe jumps in the pressure across the interface equivalent to the prescribed over-pressure forces, i.e. <math display="inline">\gamma = 5</math> for the first fluid-pair (Figure [[#img-15|15a]]) and <math display="inline">\gamma = 1000</math> for the oil&#8211;water fluid-pair (Figure [[#img-15|15c]]). Figures [[#img-15|15b]] and [[#img-15a|15d]] present the same result as shown in Figures [[#img-15|15a]] and [[#img-15|15c]], respectively, but with the solutions projected onto the <math display="inline">x-z</math> plane. In this view, all the inter-element jumps in the obtained pressure solutions are better seen.
1474
1475
<div id='img-16a'></div>
1476
<div id='img-16'></div>
1477
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1478
|-
1479
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_vol_err.png|294px|]]
1480
|[[Image:Draft_Samper_536119838-fig_ex3_umesh_tri_grid_24_12_unsteady2_rho_diff_vol_err_mod.png|294px|]]
1481
|-style="text-align: center; font-size: 75%;"
1482
|(a) Signed relative error in areas.
1483
|(b) Absolute relative error (log-scale) in areas.
1484
|-
1485
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_vol_err.png|294px|]]
1486
|[[Image:Draft_Samper_536119838-fig_ex4_umesh_tri_grid_24_12_unsteady2_rho_diff_vol_err_mod.png|294px|]]
1487
|-style="text-align: center; font-size: 75%;"
1488
|(c) Signed relative error in areas.
1489
|(d) Absolute relative error (log-scale) in areas.
1490
|- style="text-align: center; font-size: 75%;"
1491
| colspan="2" | '''Figure 16:''' Lagrangian example with fixed-walls and an unstable interface: time evolution of the area errors are shown. The top row corresponds to a fluid-pair: <math>(\rho _{1},\mu _{1}) = (1, 10^{-6})</math>, <math>(\rho _{2},\mu _{2}) = (5, 10^{-5})</math> and <math>\gamma = 5</math>. The bottom row corresponds to oil&#8211;water interface: <math>(\rho _{1},\mu _{1}) = (919, 4.9 \times 10^{-2})</math> and <math>(\rho _{2},\mu _{2}) = (997, 8.9 \times 10^{-4})</math> and <math>\gamma = 1000</math>.
1492
|}
1493
1494
The time evolution of the area-errors is shown in Figure [[#img-16|16]] (both signed and absolute relative errors). The area-errors are measured in the same way as it is done earlier in section [[#10.3 Moving-mesh (Lagrangian) extrusion problem|10.3]], i.e. using the expressions given in Eq. ([[#eq-68|68]]). For the simulated time period and the first fluid-pair, we find <math display="inline">|e^{h}_{\Omega _{1}}| < 8 \times 10^{-4}</math>, <math display="inline">|e^{h}_{\Omega _{2}}| < 8 \times 10^{-4}</math> and <math display="inline">|e^{h}_{\Omega }| < 10^{-8}</math>. Likewise for the oil&#8211;water fluid-pair we find <math display="inline">|e^{h}_{\Omega _{1}}| < 4 \times 10^{-4}</math>, <math display="inline">|e^{h}_{\Omega _{2}}| < 4 \times 10^{-4}</math> and <math display="inline">|e^{h}_{\Omega }| < 10^{-8}</math>. These results indicate the good area preservation property of the  P1-P0+ triangle.
1495
1496
==11 CONCLUDING REMARKS==
1497
1498
We have presented a 3-noded triangle and a 4-noded tetrahedra with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown ''constant pressure field'' and ''a prescribed linear field'' that satisfies the steady state momentum equations for a constant body force. In the so called P1/P0+ elements the “effective” pressure field is linear, although the unknown pressure field is piecewise constant within each element. The P1/P0+ elements  have shown an excellent behaviour in terms of accuracy and mass conservation for incompressible viscous flow problems with discontinuous material properties formulated in either Eulerian or Lagrangian descriptions.
1499
1500
==ACKNOWLEDGEMENTS==
1501
1502
This research was  supported by project HFLUIDS of the Ministerio de Educación y Ciencia  of Spain and  projects SAFECON and REALTIME of the European Research Council.
1503
1504
The authors are grateful to Prof. Shijie Zhong and Dr. Thibault Duretz for sharing their codes to evaluate the analytical solution of the buoyant Stokes flow example with a columnar viscosity structure.
1505
1506
==APPENDIX A. LINEARIZATION OF THE NAVIER-STOKES EQUATIONS==
1507
1508
Consider the steady Navier&#8211;Stokes equations described in an Eulerian framework. The stabilized momentum and  mass balance equations discretized in space by the FEM lead to the following system of equations,
1509
1510
<span id="eq-1"></span>
1511
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1512
|-
1513
| 
1514
{| style="text-align: left; margin:auto;width: 100%;" 
1515
|-
1516
| style="text-align: center;" | <math>\begin{bmatrix}(\mathbf{K}_0 + \mathbf{K}_v+\mathbf{K}_s) & \mathbf{Q}
1517
1518
\\  \mathbf{Q}^{\hbox{T}}  & \mathbf{S} \end{bmatrix} \begin{Bmatrix}\bar{\bf v}
1519
1520
\\  \bar {\bf p} \end{Bmatrix} = \begin{Bmatrix}{\bf f}_v
1521
1522
\\  {\bf f}_p \end{Bmatrix} </math>
1523
|}
1524
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
1525
|}
1526
1527
It is assumed that the stabilization terms in matrices <math display="inline">\mathbf{K}_s</math> and <math display="inline">\mathbf{S}</math> and vector <math display="inline">{\bf f}_p</math> are linear, i.e. they do not depend on the solution. The only nonlinear term that remains is, therefore,  matrix <math display="inline">\mathbf{K}_v</math> obtained from the Galerkin FEM discretization of the convective term (Eq.(55a)).
1528
1529
The Newton&#8211;Raphson linearization of Eq.([[#eq-1|1]]) gives the following system of equations.
1530
1531
<span id="eq-2"></span>
1532
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1533
|-
1534
| 
1535
{| style="text-align: left; margin:auto;width: 100%;" 
1536
|-
1537
| style="text-align: center;" | <math>\begin{bmatrix}\mathbf{H} & \mathbf{Q}
1538
1539
\\  \mathbf{Q}^{\hbox{T}}  & \mathbf{S} \end{bmatrix} \begin{Bmatrix}\delta \bar{\bf v}
1540
1541
\\  \delta \bar {\bf p} \end{Bmatrix} = \begin{Bmatrix}{\bf f}_v - (\mathbf{K}_0 + \mathbf{K}_v+\mathbf{K}_s) \bar{\bf v}^{m} - \mathbf{Q} \bar {\bf p}^{m}
1542
1543
\\  {\bf f}_p - \mathbf{Q}^{\hbox{T}} \bar{\bf v}^{m}  - \mathbf{S}_{p} \bar {\bf p}^{m} \end{Bmatrix}  </math>
1544
|}
1545
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
1546
|}
1547
1548
with
1549
1550
<span id="eq-3"></span>
1551
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1552
|-
1553
| 
1554
{| style="text-align: left; margin:auto;width: 100%;" 
1555
|-
1556
| style="text-align: center;" | <math>\mathbf{H}^e = \mathbf{K}_0^e + \mathbf{K}_v^e + \bar{\mathbf{K}}_v^e + \mathbf{K}_s^e  </math>
1557
|}
1558
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.3)
1559
|}
1560
1561
where matrices <math display="inline">\mathbf{K}_0^e,~\mathbf{K}_v^e</math> and <math display="inline">\mathbf{K}_s^e</math> are given in Eqs.([[#eq-54|54]]) and ([[#8 EXTENSION TO THE NAVIER-STOKES EQUATIONS|8]]) and
1562
1563
<span id="eq-4"></span>
1564
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1565
|-
1566
| 
1567
{| style="text-align: left; margin:auto;width: 100%;" 
1568
|-
1569
| style="text-align: center;" | <math>\bar {\bf K}^e_{v_{ij}} =  \int _{\Omega ^e}\rho \mathbf{N}_i  ({\boldsymbol \nabla } {v}^T)^T \mathbf{N}_jd\Omega   </math>
1570
|}
1571
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.4)
1572
|}
1573
1574
In Eq.([[#eq-2|2]]) index <math display="inline">m</math>  denotes the iterations. Hence,
1575
1576
<span id="eq-5"></span>
1577
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1578
|-
1579
| 
1580
{| style="text-align: left; margin:auto;width: 100%;" 
1581
|-
1582
| style="text-align: center;" | <math>\bar{\bf v}^{m+1} =\bar{\bf v}^{m} +\delta \bar{\bf v}\quad \hbox{and}\quad \bar{\bf p}^{m+1} =\bar{\bf p}^{m} +\delta \bar{\bf p}  </math>
1583
|}
1584
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.5)
1585
|}
1586
1587
The Picard linearization of Eq.([[#eq-1|1]]) yields the following system of equations
1588
1589
<span id="eq-6"></span>
1590
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1591
|-
1592
| 
1593
{| style="text-align: left; margin:auto;width: 100%;" 
1594
|-
1595
| style="text-align: center;" | <math>\begin{bmatrix}(\mathbf{K} + \mathbf{K}_v^m + \mathbf{K}_s) & \mathbf{Q}
1596
1597
\\  \mathbf{Q}^{\hbox{T}} & \mathbf{S} \end{bmatrix} \begin{Bmatrix}\bar{\bf v}^{m+1}
1598
1599
\\  \bar {\bf p}^{m+1} \end{Bmatrix} = \begin{Bmatrix}{\bf f}_v
1600
1601
\\  {\bf f}_p \end{Bmatrix} </math>
1602
|}
1603
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.6)
1604
|}
1605
1606
==APPENDIX B. TIME INTEGRATION IN THE LAGRANGIAN FRAMEWORK==
1607
1608
The momentum and the stabilized mass balance equations discretized in space by the FEM leads to the following system of equations in the Lagrangian framework
1609
1610
<span id="eq-1"></span>
1611
<span id="eq-2"></span>
1612
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1613
|-
1614
| 
1615
{| style="text-align: left; margin:auto;width: 100%;" 
1616
|-
1617
| style="text-align: center;" | <math>\mathbf{M}\bar{\bf v} + \mathbf{K} \bar{\bf v} + \mathbf{Q} \bar {\bf p} = {\bf f}_v </math>
1618
| style="width: 5px;text-align: right;white-space: nowrap;" | (B.1)
1619
|-
1620
| style="text-align: center;" | <math>  \mathbf{Q}^{\hbox{T}}  \bar{\bf v} + \mathbf{S} \bar {\bf p} = {\bf f}_p </math>
1621
| style="width: 5px;text-align: right;white-space: nowrap;" | (B.2)
1622
|}
1623
|}
1624
1625
where, the matrices and vectors with sub-scripts are the ones that appear due to the stabilization of the mass balance equation. Let <math display="inline">\bar{a}= {\partial \bar v \over \partial t}</math> denote the vector of unknown nodal accelerations. The momentum balance equation, i.e. Eq.([[#eq-1|B.1]]), is a dynamic equation. On the contrary, the stabilized mass balance equation, i.e. Eq.([[#eq-2|B.2]]), is a quasi-static equation.
1626
1627
Using the Newmark algorithm the balance equations given in Eqs.([[#eq-1|B.1]]) and ([[#eq-2|B.2]]) are integrated in time (<math display="inline">\{ {}^{n}\mathbf{x}, {}^{n}\mathbf{v}, {}^{n}\mathbf{a}, {}^{n}\bar {p}\} \rightarrow \{ {}^{n+1} {x}, {}^{n+1}\bar{v}, {}^{n+1}\bar{a}, {}^{n+1}\bar {p}\} </math>) as follows
1628
1629
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1630
|-
1631
| 
1632
{| style="text-align: left; margin:auto;width: 100%;" 
1633
|-
1634
| style="text-align: center;" | <math>\begin{bmatrix}\dfrac{1}{\Delta t} \mathbf{M} + \theta ~{}^{n+1}\mathbf{K} & \theta ~{}^{n+1}\mathbf{Q}
1635
1636
\\  \theta \mathbf{Q}^{\hbox{T}}  & \theta ~{}^{n+1}\mathbf{S} \end{bmatrix} \begin{Bmatrix}{}^{n+1}\bar{\bf v}
1637
1638
\\  {}^{n+1}\bar {\bf p} \end{Bmatrix} = \begin{Bmatrix}\theta ~{}^{n+1} {\bf f}_v + {}^{n}{\bf r}_m
1639
1640
\\  \theta  ~{}^{n+1} {\bf f}_p \end{Bmatrix} </math>
1641
|-
1642
| style="text-align: center;" | <math> {}^{n}{\bf r}_m := \frac{1}{\Delta t} \mathbf{M} {}^{n}\bar{\bf v}- (1-\theta ) \left[{}^{n}\mathbf{K} {}^{n}\bar{\bf v} + {}^{n}\mathbf{Q} {}^{n}\bar {\bf p} - {}^{n}{\bf f}_v \right]</math>
1643
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
1644
|-
1645
| style="text-align: center;" | <math> {}^{n+1}\bar{\bf a} = \frac{1}{\theta \Delta t} \left(\bar{\bf v}^{n+1} - {}^{n}\bar{\bf v} \right)- \frac{(1-\theta )}{\theta } {}^{n}{\bf a} </math>
1646
|-
1647
| style="text-align: center;" | <math> {}^{n+1}{\bf x} = {}^{n}{\bf x} + {\Delta t} {}^{n}\bar{\bf v} + \frac{\Delta t^{2}}{2} \left[(1-2\beta ){}^{n}\bar{\bf a} + 2\beta {}^{n+1}\bar{\bf a} \right] </math>
1648
|}
1649
|}
1650
1651
==References==
1652
1653
<div id="cite-1"></div>
1654
'''[[#citeF-1|[1]]]'''  R. Aubry, S.R. Idelsohn, E. Oñate,  Particle finite element method in fluid mechanics including thermal convection&#8211;diffusion, Comput. Struct. 83 (2004) 1459–-75.
1655
1656
<div id="cite-2"></div>
1657
'''[[#citeF-2|[2]]]'''   R.F. Ausas , G.C. Buscaglia, S.R. Idelsohn, A new enrichment space for the treatment of discontinuous pressures in multi-fluid flows, Int. J. Numer. Meth. Fluids 70 (7) (2012) 829&#8211;850.
1658
1659
<div id="cite-3"></div>
1660
'''[[#citeF-3|[3]]]'''  R. Araya, G.R. Barrechea, F. Valentin,  Stabilized finite element methods based on multiscale enrichment for the Stokes problem, Siam J. Numer. Anal. 44 (1) (2006) 322–-348.
1661
1662
<div id="cite-4"></div>
1663
'''[[#citeF-4|[4]]]'''  T. Belytschko, W.K. Liu, B. Moran, Non linear finite element for continua and structures, 2d Edition, Wiley, 2013.
1664
1665
<div id="cite-5"></div>
1666
'''[[#citeF-5|[5]]]'''  E. Burman, M. Fernández, P. Hansbo,  Edge stabilizatio n for the incompressible Navier–Stokes equations: a continuous interior penalty finite element method, Tech. Report RR-5349, INRIA, Le Chesnay, France, 2004.
1667
1668
<div id="cite-6"></div>
1669
'''[[#citeF-6|[6]]]'''   E. Burman,  P. Hansbo,  A unified stabilized method for Stokes’ and Darcy equations, Tech. Report 2002-15, Chalmers Finite Element Center, Göteborg, Sweden, 2002.
1670
1671
<div id="cite-7"></div>
1672
'''[[#citeF-7|[7]]]'''  A. Caiazzo, M.A. Fernández, V. Martin, Analysis of a stabilized finite element method for fluid flows through a porous interface, Applied Mathematics Letters 24 (2011) 2124&#8211;2127.
1673
1674
<div id="cite-8"></div>
1675
'''[[#citeF-8|[8]]]'''  J.M. Carbonell, E. Oñate, B. Suárez, 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.
1676
1677
<div id="cite-9"></div>
1678
'''[[#citeF-9|[9]]]'''  R. Codina, S. Badia,  On the design of discontinuous Galerkin methods for elliptic problems based on hybrid formulation, Comput. Methods Appl. Mech. Engrg. (2013) http://dx.doi.org/10.1016/j.cma.2013.05.004.
1679
1680
<div id="cite-10"></div>
1681
'''[[#citeF-10|[10]]]'''  F. Del Pin, S.R. Idelsohn, E. Oñate, R. Aubry, The ALE/Lagrangian Particle Finite Element Method: A new approach to computation of free-surface flows and fluid–object interactions, Computers & Fluids  36 (1) (2007) 27&#8211;38.
1682
1683
<div id="cite-11"></div>
1684
'''[[#citeF-11|[11]]]'''   J. Donea,  A. Huerta,   Finite element method for flow problems, J. Wiley & Sons, 2003, pp. 362.
1685
1686
<div id="cite-12"></div>
1687
'''[[#citeF-12|[12]]]'''  J. Douglas, J. Wang, An absolutely stabilized finite element method for the Stokes problem,  Math. Comp. 52 (186) (1989) 495&#8211;508.
1688
1689
<div id="cite-13"></div>
1690
'''[[#citeF-13|[13]]]'''  T.&nbsp;Duretz, D.&nbsp;A. May, T.&nbsp;V. Gerya, P.&nbsp;J. Tackley, Discretization errors and   free surface stabilization in the finite difference and marker-in-cell method   for applied geodynamics: A numerical study, Geochemistry, Geophysics,   Geosystems 12&nbsp;(7) (2011).
1691
1692
<div id="cite-14"></div>
1693
'''[[#citeF-14|[14]]]'''  M.A. Fernández, J.-F. Gerbeau, V. Martin,  Numerical simulation of blood flows through a porous interface, ESAIM: Mathematical Modelling and Numerical Analysis 42 (2008) 961&#8211;990.
1694
1695
<div id="cite-15"></div>
1696
'''[[#citeF-15|[15]]]'''  F.&nbsp;Gilbert, G.&nbsp;E. Backus, Propagator matrices in elastic wave and vibration problems, Geophysics 31&nbsp;(2) (1966) 326&#8211;332.
1697
1698
<div id="cite-16"></div>
1699
'''[[#citeF-16|[16]]]'''  B.&nbsp;H. Hager, R.&nbsp;J. O'Connell, A simple global model of   plate dynamics and mantle convection, Journal of Geophysical Research   86&nbsp;(B6) (1981) 4843&#8211;4867.
1700
1701
<div id="cite-17"></div>
1702
'''[[#citeF-17|[17]]]'''  T.J.R. Hughes,  L.P. Franca,   A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 85–-96.
1703
1704
<div id="cite-18"></div>
1705
'''[[#citeF-18|[18]]]'''  S.R. Idelsohn,  E. Oñate,  F. Del Pin,   The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves, Int. J. Numer. Methods Engrg.  61 (2004) 964&#8211;989.
1706
1707
<div id="cite-19"></div>
1708
'''[19]'''  S.R. Idelsohn,  E. Oñate,  F. Del Pin,  Fluid-structure interaction with the Particle Finite Element Method, Comput. Methods Appl. Mech. Engrg. 195 (2006) 2100&#8211;2113.
1709
1710
<div id="cite-20"></div>
1711
'''[[#citeF-20|[20]]]'''  S.R. Idelsohn, M. Mier-Torrecilla, E. Oñate,  Multi-fluid flows with the Particle Finite Element Method, Comput Methods Appl Mech Engrg. 198  (2009) 2750&#8211;2767.
1712
1713
<div id="cite-21"></div>
1714
'''[[#citeF-21|[21]]]'''  S.R. Idelsohn, M. Mier-Torrecilla, R. Nigro, E. Oñate,   On the analysis of heterogeneous fluids with jumps in the viscosity using a discontinuous pressure field, Comput. Mech.  46 (1) (2010) 115&#8211;124
1715
1716
<div id="cite-22"></div>
1717
'''[[#citeF-22|[22]]]'''  S.R. Idelsohn,  E. Oñate,  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 (2010) 1313–-1330.
1718
1719
<div id="cite-23"></div>
1720
'''[[#citeF-23|[23]]]'''  N. Kechar, D. Silvester,  Analysis of a locally stabilized mixed finite element method for the Stokes problem, Math. Comp. 58 (1992) 1–-10.
1721
1722
<div id="cite-24"></div>
1723
'''[[#citeF-24|[24]]]'''  M.&nbsp;Kronbichler, T.&nbsp;Heister, W.&nbsp;Bangerth, High   accuracy mantle convection simulation through modern numerical methods,   Geophysical Journal International 191&nbsp;(1) (2012) 12&#8211;29.
1724
1725
<div id="cite-25"></div>
1726
'''[[#citeF-25|[25]]]'''  M. Mier-Torrecilla, S.R. Idelsohn, E. Oñate, Advances in the simulation of multi-fluid flows with the particle finite element method. Application to bubble dynamics, Int. J. Numer. Meth. Fluids 67 (2011) 1516–-1539.
1727
1728
<div id="cite-26"></div>
1729
'''[[#citeF-26|[26]]]'''  E. Oñate,   Derivation of stabilized equations for advective-diffusive transport and fluid flow problems, Comput. Meth. Appl. Mech. Engng. 151 (1998) 233&#8211;267.
1730
1731
<div id="cite-27"></div>
1732
'''[27]'''   E. Oñate, A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation, Comput Methods Appl Mech Engrg. 182 (1&#8211;2) (2000) 355&#8211;370.
1733
1734
<div id="cite-28"></div>
1735
'''[[#citeF-28|[28]]]'''   E. Oñate, J. García,   A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation, Comput. Meth. Appl. Mech. Engrg. 191 (2001) 635&#8211;660.
1736
1737
<div id="cite-29"></div>
1738
'''[29]'''  E. Oñate,  Multiscale computational analysis in mechanics using   finite calculus: an introduction, Comput. Meth. Appl. Mech. Engrg.   192(28-30) (2003) 3043&#8211;3059.
1739
1740
<div id="cite-30"></div>
1741
'''[[#citeF-30|[30]]]'''   E. Oñate,  Possibilities of finite calculus in computational mechanics, Int. J. Num. Meth. Engng. 60 (1) (2004) 255&#8211;281.
1742
1743
<div id="cite-31"></div>
1744
'''[[#citeF-31|[31]]]'''   E. Oñate, J. Rojek, R.L. Taylor, O.C. Zienkiewicz,  Finite calculus formulation for incompressible solids using linear triangles and tetrahedra, Int. J. Numer. Meth. Engng. 59 (11) (2004a) 1473&#8211;1500.
1745
1746
<div id="cite-32"></div>
1747
'''[[#citeF-32|[32]]]'''  E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry, The particle finite element method. An overview, Int. J. Comput. Methods 1 (2) (2004b) 267&#8211;307.
1748
1749
<div id="cite-33"></div>
1750
'''[[#citeF-33|[33]]]'''  E. Oñate, A. Valls, J. García,   FIC/FEM formulation   with matrix stabilizing terms for incompressible flows at low and high   Reynold's numbers, Comput. Mech. 38 (4-5) (2006) 440&#8211;455.
1751
1752
<div id="cite-34"></div>
1753
'''[34]'''  E. Oñate, J. García, S.R. Idelsohn, F. Del Pin, Finite calculus formulation for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches, Comput. Meth. Appl. Mech. Engng. 195 (2006b) 3001&#8211;3037.
1754
1755
<div id="cite-35"></div>
1756
'''[[#citeF-35|[35]]]'''  E. Oñate, A. Valls, J. García,   Computation of turbulent flows using a finite calculus-finite element formulation, Int. J. Numer. Meth. Engng. 54 (2007) 609&#8211;637.
1757
1758
<div id="cite-36"></div>
1759
'''[[#citeF-36|[36]]]'''  E. Oñate, M.A. Celigueta, S.R. Idelsohn,  F. Salazar, B. Suárez, Possibilities of the particle finite element method for fluid-soil-structure interaction problems, Comput. Mech. 48(3) (2011)  307&#8211;318.
1760
1761
<div id="cite-37"></div>
1762
'''[[#citeF-37|[37]]]'''  E. Oñate, A. Franci, J.M. Carbonell, A FIC-based stabilized Lagrangian formulation with negligible mass losses for incompressible fluids. Application to free-surface flows using PFEM, Publication CIMNE, No. PI394. Submitted to Int. J. Num. Meth. Fluids, 2013.
1763
1764
<div id="cite-38"></div>
1765
'''[[#citeF-38|[38]]]'''  P. Ryzhakov, R. Rossi, S.R. Idelsohn, E. Oñate, A monolithic Lagrangian approach for fluid-structure interaction problems, Comput. Mech. 46 (6) 2010 883&#8211;899.
1766
1767
<div id="cite-39"></div>
1768
'''[[#citeF-39|[39]]]'''  O.C. Zienkiewicz,  R.L. Taylor, P. Nithiarasu,  The finite element   method for fluid dynamics,  Elsevier, 2006.
1769
1770
<div id="cite-40"></div>
1771
'''[[#citeF-40|[40]]]'''  S.&nbsp;Zhong. Analytic   solutions for Stokes’ flow with lateral variations in viscosity, ''Geophysical Journal International'' 124&nbsp;(1) 1996 18&#8211;28.
1772

Return to Onate et al 2015a.

Back to Top

Document information

Published on 01/01/2015

DOI: 10.1016/j.cma.2013.12.009
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 4
Views 26
Recommendations 0

Share this document