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
<!-- metadata commented in wiki content
2
==STABILIZATION VIA FINITE CALCULUS IN FINITEELEMENT COMPUTATIONS==
3
4
''' 
5
Eugenio Oñate  '''
6
7
{|
8
|-
9
|   International Center for Numerical 
10
|-
11
| Methods in Engineering (CIMNE) 
12
|-
13
| Universidad Politecnica de Cataluna 
14
|-
15
| Campus Norte UPC, 08034 Barcelona, Spain 
16
|-
17
| Email: [mailto:onate@cimne.upc.es onate@cimne.upc.es] 
18
|-
19
| Web page: [http://www.cimne.upc.es http://www.cimne.upc.es] 
20
|-
21
| Dedicated to Tom Hughes on the occasion of his 60th birthday   
22
|}
23
-->
24
25
==Abstract==
26
27
The expression “finite calculus” refers to the derivation of the governing differential equations in mechanics by invoking balance of fluxes, forces, etc. in a space-time domain of finite size. The governing equations resulting from this approach are different from those of infinitesimal calculus theory and they incorporate new terms which depend on the dimensions of the balance domain. The new governing equations allow to derive naturally stabilized numerical schemes using any discretization procedure. The paper  discusses the possibilities of the  finite calculus method  for the finite element solution of convection-diffusion problems with sharp gradients and incompressible fluid flow.
28
29
'''Keywords''' Stabilization, finite calculus, finite element method.
30
31
==1 INTRODUCTION==
32
33
It is well known that standard numerical methods such as the central finite difference (FD) method, the Galerkin finite element (FE) method and the finite volume (FV) method, among others, lead to unstable numerical solutions when applied to problems involving different scales, multiple constraints and/or high gradients. Examples of these situations are typical in the solution of convection-diffusion problems, incompressible problems in fluid and solid mechanics and strain or strain rate localization problems in solids and compressible fluids using the standard Galerkin FE method or central scheme in FD and  FV methods <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]]]. Similar instabilities are found in the application of meshless methods to those problems <span id='citeF-3'></span> <span id='citeF-4'></span> <span id='citeF-5'></span>[[#cite-3|[3]]-[[#cite-5|5]]].
34
35
The sources of the numerical instabilities in FE, FD and FV methods, for instance, have been sought in the apparent unability of the Galerkin FE method and the analogous central difference scheme in FD and FV methods, to provide a numerical procedure able to capture the different scales appearing in the  solution for all ranges of the physical parameters. Typical examples  are the spurious numerical oscillations  in convection-diffusion problems for high values of the convective terms. The same type of  oscillations are found in regions next to  sharp internal layers appearing in high speed compressible flows (shocks) or in strain localization problems (shear bands) in solids. A similar problem of different nature emerges in the solution of incompressible  problems in fluid and solid mechanics. Here the difficulties in satisfying the incompressibility constraint limit the choices of the approximation for the velocity (or displacement) variables and the pressure <span id='citeF-1'></span>[[#cite-1|[1]]].
36
37
The solution of above problems has been attempted in a number of ways. The underdiffusive character of the central difference scheme for treating advective-diffusive problems has been corrected in an ad-hoc manner by adding the so called “artificial diffusion” terms to the standard governing equation <span id='citeF-2'></span>[[#cite-2|[2]]]. The same idea has been successfully applied to derive stabilized FV and FE methods for convection-diffusion and fluid-flow problems <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]]]. Other stabilized FD schemes are based on the “upwind” computation of the first derivatives appearing in the convective operator <span id='citeF-2'></span>[[#cite-2|[2]]]. The counterpart of upwind techniques in the FEM are the so called Petrov-Galerkin methods <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-6'></span>[[#cite-6|6]],<span id='citeF-7'></span>[[#cite-7|7]]]. Among the many methods of this kind we can name the SUPG method <span id='citeF-8'></span><span id='citeF-9'></span>[[#cite-8|[8]]-<span id='citeF-10'></span>[[#cite-10|10]]], the Galerkin Least Square (GLS) method <span id='citeF-11'></span>[[#cite-11|[11]],<span id='citeF-12'></span>[[#cite-12|12]]] the Characteristic Galerkin method <span id='citeF-13'></span><span id='citeF-14'></span>[[#cite-13|[13]]-<span id='citeF-15'></span>[[#cite-15|15]]] the Characteristic Based Split (CBS) method <span id='citeF-16'></span>[[#cite-16|[16]],<span id='citeF-17'></span>[[#cite-17|17]]] and the Subgrid Scale (SS) method <span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-18|[18]]-<span id='citeF-21'></span>[[#cite-21|21]]].
38
39
The Finite Calculus (FIC) is a  different route to derive stabilized numerical methods. The starting point are the modified governing differential equations of the problem derived by expressing the balance of fluxes (or equilibrium of forces) in a space-time domain of finite size <span id='citeF-22'></span>[[#cite-22|[22]]]. This introduces naturally new terms in the classical differential equations of the infinitesimal theory which are a function of the balance domain dimensions. The merit of the modified equations via the FIC approach is that they lead to stabilized schemes using ''any'' numerical method. In addition, the different stabilized FD, FE and FV methods typically used in practice can be ''recovered'' using the FIC  equations <span id='citeF-22'></span>[[#cite-22|[22]],<span id='citeF-23'></span>[[#cite-23|23]]]. Moreover, these equations are the basis for deriving a  procedure for computing the stabilization parameters <span id='citeF-24'></span>[[#cite-24|[24]],<span id='citeF-25'></span>[[#cite-25|25]]].
40
41
Most stabilized FEM schemes can be framed within an extended Galerkin approach where the standard Galerkin expression is modified by adding adequate residual-based terms as
42
43
<span id='eq-1'></span>
44
{| class="formulaSCP" style="width: 100%; text-align: left;" 
45
|-
46
| 
47
{| style="text-align: left; margin:auto;width: 100%;" 
48
|-
49
| style="text-align: center;" | <math>[\hbox{Galerkin}] +\sum \limits _e \int _{\Omega ^e} [\tau ^e] {\boldsymbol P}  (N_k) \cdot {\boldsymbol r} d\Omega =0 </math>
50
|}
51
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
52
|}
53
54
where [Galerkin] denotes the standard Galerkin expression, <math display="inline">{\boldsymbol P}(N_j)</math> is a vector which terms depend on the shape functions <math display="inline">N_j</math> (and the physical parameters of the problem), '''r''' is the vector of residuals of the finite element equations and <math display="inline">[\tau ^e]</math> is a matrix of stabilization parameters.
55
56
The computation of the stabilization parameters is still an open problem and much effort has been devoted to this topic <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]],<span id='citeF-7'></span>[[#cite-7|7]]-<span id='citeF-21'></span>[[#cite-21|21]],<span id='citeF-36'></span>[[#cite-36|36]]].
57
58
Different stabilized methods can be derived from Eq.([[#eq-1|1]]) by chosing expressions for matrices '''P''' and <math display="inline">[\tau ^e]</math>.  In this paper we will focus in the stabilized finite element formulation using the FIC method.  The equivalent form of Eq.([[#eq-1|1]]) is written in the FIC formulation as <span id='citeF-1'></span>[[#cite-1|[1]]-<span id='citeF-7'></span>[[#cite-7|7]]]
59
60
<span id='eq-2'></span>
61
{| class="formulaSCP" style="width: 100%; text-align: left;" 
62
|-
63
| 
64
{| style="text-align: left; margin:auto;width: 100%;" 
65
|-
66
| style="text-align: center;" | <math>\int _\Omega N_k r_i d\Omega + \sum \limits _e \int _{\Omega ^e} {h_j\over 2} {\partial N_k\over \partial x_j} r_i d\Omega =0\quad i=1,n_r </math>
67
|}
68
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
69
|}
70
71
where <math display="inline">n_{r}</math> is the number of residual equations, <math display="inline">r_i</math> is the <math display="inline">i</math>th residual equation and <math display="inline">h_j</math> are characteristic length parameters which are typically of the  order of the element dimensions.
72
73
The FIC formulation has been used in conjunction with the finite element formulation to solve a variety of problems in convection-diffusion <span id='citeF-23'></span>[[#cite-23|[23]]-<span id='citeF-27'></span>[[#cite-27|27]]] incompressible fluid dynamics involving free surfaces <span id='citeF-28'></span>[[#cite-28|[28]]-<span id='citeF-32'></span>[[#cite-32|32]]] and non linear solid mechanics problems allowing for large strains <span id='citeF-28'></span>[[#cite-28|[28]],<span id='citeF-41'></span>[[#cite-41|41]],<span id='citeF-42'></span>[[#cite-42|42]]] using in all cases linear triangles and tetrahedra with equal interpolation for all variables.
74
75
The layout of the paper is the following. In the next section the main concepts of the FIC method are introduced. Applications of the FIC method to convection-diffusion problems with sharp gradients are detailed and some examples of application are given. Finally the possibilities of the FIC method in incompressible fluid mechanics are discussed and a finite element formulation is presented.
76
77
==2 THE FINITE CALCULUS METHOD==
78
79
We will consider a convection-diffusion problem in a 1D domain <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure [[#img-1|1]]) is written as
80
81
<span id='eq-3'></span>
82
{| class="formulaSCP" style="width: 100%; text-align: left;" 
83
|-
84
| 
85
{| style="text-align: left; margin:auto;width: 100%;" 
86
|-
87
| style="text-align: center;" | <math>q_A - q_B=0 </math>
88
|}
89
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
90
|}
91
92
where <math display="inline">q_A</math> and <math display="inline">q_B</math> are the incoming and outgoing fluxes at points <math display="inline">A</math> and <math display="inline">B</math>, respectively. The flux <math display="inline">q</math> includes both convective and diffusive terms; i.e. <math display="inline">q=v\phi - k{d\phi \over dx}</math>, where <math display="inline">\phi </math> is the transported variable (i.e. the temperature in a thermal problem), <math display="inline">v</math> is the velocity and <math display="inline">k</math> is the diffusitivity of the material.
93
94
<div id='img-1'></div>
95
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
96
|-
97
|[[Image:Draft_Samper_674218428-majesus2.png|300px|Equilibrium of fluxes in a  space balance domain of finite size]]
98
|- style="text-align: center; font-size: 75%;"
99
| colspan="1" | '''Figure 1:''' Equilibrium of fluxes in a  space balance domain of finite size
100
|}
101
102
Let us express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure [[#img-1|1]]). Expanding <math display="inline">q_A</math> and <math display="inline">q_B</math> in Taylor series around point <math display="inline">C</math> up to second order terms gives
103
104
<span id='eq-4'></span>
105
{| class="formulaSCP" style="width: 100%; text-align: left;" 
106
|-
107
| 
108
{| style="text-align: left; margin:auto;width: 100%;" 
109
|-
110
| style="text-align: center;" | <math>q_A= q_C - d_1 \frac{dq}{d x}\vert _C+ \frac{d^2_1}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_1)\quad ,\quad  q_B= q_C + d_2 \frac{dq}{d x}\vert _C+\frac{d^2_2}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_2) </math>
111
|}
112
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
113
|}
114
115
Substituting Equation ([[#eq-2|2]]) into Equation ([[#eq-1|1]]) gives after simplification
116
117
<span id='eq-5'></span>
118
{| class="formulaSCP" style="width: 100%; text-align: left;" 
119
|-
120
| 
121
{| style="text-align: left; margin:auto;width: 100%;" 
122
|-
123
| style="text-align: center;" | <math>\frac{dq}{dx}-\underline{\frac{h}{2} \frac{d^2q}{dx^2}}=0 </math>
124
|}
125
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
126
|}
127
128
where <math display="inline">h=d_1-d_2</math> and all the derivatives are computed at the arbitrary point <math display="inline">C</math>.
129
130
Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the space balance domain to have a ''finite size''. The new balance equation ([[#eq-3|3]]) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in Equation ([[#eq-2|2]]) would lead to new terms in Equation ([[#eq-3|3]]) involving higher powers of <math display="inline">h</math>.
131
132
Distance <math display="inline">h</math> in Equation ([[#eq-3|3]]) can be interpreted as a free parameter depending on the location of point <math display="inline">C</math> within the balance domain. Note that <math display="inline">-d\le h\le d</math> and, hence, <math display="inline">h</math> can take a negative value. At the discrete solution level the domain <math display="inline">d</math> should be replaced by the balance domain around a node. This gives for an equal size discretization <math display="inline">-l^e\le h \le l^e</math> where <math display="inline">l^e</math> is the element or cell dimension. The fact that Equation ([[#eq-3|3]]) is the ''exact balance equation'' (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule leading to an smaller error in the numerical solution <span id='citeF-22'></span>[[#cite-22|[22]]-<span id='citeF-24'></span>[[#cite-24|24]]].
133
134
Consider, for instance, Equation ([[#eq-3|3]]) applied to the 1D convection-diffusion problem. Neglecting third order derivatives of <math display="inline">\phi </math>, Equation ([[#eq-3|3]]) can be rewritten  in terms of <math display="inline">\phi </math> as
135
136
<span id='eq-6'></span>
137
{| class="formulaSCP" style="width: 100%; text-align: left;" 
138
|-
139
| 
140
{| style="text-align: left; margin:auto;width: 100%;" 
141
|-
142
| style="text-align: center;" | <math>-v \frac{d\phi }{dx}+\left(k+\frac{v h}{2}\right)\frac{d^2\phi }{dx^2}=0 </math>
143
|}
144
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
145
|}
146
147
We see clearly that the FIC method introduces ''naturally'' an additional diffusion term in the standard convection-diffusion equation. This is the basis of the popular “artificial diffusion” procedure <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]],<span id='citeF-6'></span>[[#cite-6|6]]] where the characteristic length <math display="inline">h</math> is typically expressed as a function of the cell or element dimension. The ''critical''  value of <math display="inline">h</math> can be computed by introducing an optimality condition, such as obtaining a physically meaningful solution. An interpretation of the FIC equations as a modified residual method is presented in <span id='citeF-28'></span>[[#cite-28|[28]]].
148
149
Equation ([[#eq-3|3]]) can be extended to account for source terms. The modified governing equation can then be  written in compact form as
150
151
<span id='eq-7'></span>
152
{| class="formulaSCP" style="width: 100%; text-align: left;" 
153
|-
154
| 
155
{| style="text-align: left; margin:auto;width: 100%;" 
156
|-
157
| style="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}=0 </math>
158
|}
159
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
160
|}
161
162
with
163
164
<span id='eq-8'></span>
165
{| class="formulaSCP" style="width: 100%; text-align: left;" 
166
|-
167
| 
168
{| style="text-align: left; margin:auto;width: 100%;" 
169
|-
170
| style="text-align: center;" | <math>r : = -v \frac{d\phi }{dx}+ \frac{d}{dx}\left(k \frac{d\phi }{dx}\right)+ Q </math>
171
|}
172
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
173
|}
174
175
where <math display="inline">Q</math> is the external source.
176
177
The essential (Dirichlet) boundary condition for  eqn.&nbsp;([[#eq-7|7]])  is the standard one (i.e. <math display="inline">\phi =\bar \phi </math> on <math display="inline">\Gamma _\phi </math> where <math display="inline">\Gamma _\phi </math> is the boundary whete the prescribed value <math display="inline">\bar \phi </math> is imposed).  For consistency a stabilized Neumann boundary condition must be obtained.
178
179
<div id='img-2'></div>
180
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
181
|-
182
|[[Image:Draft_Samper_674218428-Neumann.png|300px|Balance domain next to a Neumann boundary point B]]
183
|- style="text-align: center; font-size: 75%;"
184
| colspan="1" | '''Figure 2:''' Balance domain next to a Neumann boundary point B
185
|}
186
187
The length of the balance segment <math display="inline">AB</math> next to a Neumann boundary is taken  as one half of the characteristic length <math display="inline">h</math> for the interior domain (Figure [[#img-2|2]]). The balance equation, assuming a constant  distribution for the source <math display="inline">Q</math> over <math display="inline">AB</math>, is
188
189
<span id='eq-9'></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>\bar q - q(x_A) - [u \phi ]_{A} - \frac{h}{2} Q = 0  </math>
196
|}
197
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
198
|}
199
200
where <math display="inline">\bar q</math> is the prescribed total flux at <math display="inline">x=L</math> and <math display="inline">x_A=x_B-\frac{h}{2}</math>.
201
202
Using a second order expansion for the advective and diffusive fluxes at point <math display="inline">A</math> gives <span id='citeF-22'></span>[[#cite-22|[22]]]
203
204
<span id="eq-10"></span>
205
{| class="formulaSCP" style="width: 100%; text-align: left;" 
206
|-
207
| 
208
{| style="text-align: left; margin:auto;width: 100%;" 
209
|-
210
| style="text-align: center;" | <math>-u \phi + k \frac{\mathrm{d} \phi }{\mathrm{d} x} + \bar q - \underline{\frac{h}{2} r} \qquad on \quad x=L </math>
211
|}
212
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
213
|}
214
215
where <math display="inline">r</math> is given by  eqn.&nbsp;([[#eq-8|8]]).
216
217
Note that for <math display="inline">h=0</math> the infinitesimal form of the 1D Neumann boundary condition is obtained.
218
219
The underlined terms in Equations ([[#eq-7|7]]) and ([[#eq-10|10]]) introduce the necessary stabilization in the discrete solution  using whatever numerical scheme.
220
221
The time dimension can be simply accounted for the FIC method by considering the balance equation in a space-time slab domain. Application of the FIC method to the transient convection-diffusion equations and to fluid flow problems can be found in <span id='citeF-23'></span>[[#cite-23|[23]],<span id='citeF-26'></span>[[#cite-26|26]],<span id='citeF-29'></span>[[#cite-29|29]]-<span id='citeF-33'></span>[[#cite-33|33]]]. Quite generally the FIC equation can be written for any problem in mechanics as <span id='citeF-22'></span>[[#cite-22|[22]]]
222
223
<span id='eq-11'></span>
224
{| class="formulaSCP" style="width: 100%; text-align: left;" 
225
|-
226
| 
227
{| style="text-align: left; margin:auto;width: 100%;" 
228
|-
229
| style="text-align: center;" | <math>r_i - \underline{\frac{h_j}{2}{\partial r_i \over \partial x_j}}-\underline{\frac{\delta }{2}{\partial r_i \over \partial t}}=0 \quad , \begin{array}{l}i=1,n_b\\ j=1,n_d \end{array} </math>
230
|}
231
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
232
|}
233
234
where <math display="inline">r_i</math> is the ith standard differential equation of the infinitesimal theory, <math display="inline">h_j</math> are characteristic length parameters, <math display="inline">\delta </math> is a time stabilization parameter and <math display="inline">t</math> the time; <math display="inline">n_b</math> and <math display="inline">n_d</math> are respectively the number of balance equations and the number of space dimensions of the problem (i.e., <math display="inline">n_d =2</math> for 2D problems, etc.). Indeed for the transient case the initial boundary conditions must be specified. The usual sum convention for repeated indexes is used in the text unless otherwise specified.
235
236
For example, in the case of the convection-diffusion problem <math display="inline">n_b=1</math> and Equation ([[#eq-8|8]]) is particularized as
237
238
<span id='eq-12'></span>
239
{| class="formulaSCP" style="width: 100%; text-align: left;" 
240
|-
241
| 
242
{| style="text-align: left; margin:auto;width: 100%;" 
243
|-
244
| style="text-align: center;" | <math>r- \underline{\frac{h_j}{2}{\partial r \over \partial x_j}}-\underline{\frac{\delta }{2}{\partial r \over \partial t}}=0\quad , \qquad j=1,n_d </math>
245
|}
246
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
247
|}
248
249
with
250
251
<span id='eq-13'></span>
252
{| class="formulaSCP" style="width: 100%; text-align: left;" 
253
|-
254
| 
255
{| style="text-align: left; margin:auto;width: 100%;" 
256
|-
257
| style="text-align: center;" | <math>r:=-\left(\displaystyle {\partial \phi  \over \partial t}+v_j \displaystyle {\partial \phi  \over \partial x_j}\right)+ \displaystyle \frac{d}{dx_j}\left(k \displaystyle \frac{d\phi }{dx_j}\right)+ Q </math>
258
|}
259
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
260
|}
261
262
For a transient solid mechanics problems Equation ([[#eq-11|11]]) applies with <math display="inline">n_b=n_d</math> and
263
264
<span id='eq-14'></span>
265
{| class="formulaSCP" style="width: 100%; text-align: left;" 
266
|-
267
| 
268
{| style="text-align: left; margin:auto;width: 100%;" 
269
|-
270
| style="text-align: center;" | <math>r_i:= -\rho \frac{\partial ^2u_i}{\partial t^2}+ {\partial \sigma _{ij} \over \partial x_j}+b_i\quad , \quad i,j=1,n_d </math>
271
|}
272
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
273
|}
274
275
where <math display="inline">u_i</math> are the displacements, <math display="inline">\sigma _{ij}</math> are the stresses and <math display="inline">b_i</math> the external body forces.
276
277
The modified Neumann boundary conditions in the FIC formulation can be written in the general case as <span id='citeF-22'></span>[[#cite-22|[22]]]
278
279
<span id='eq-15'></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>q_{ij}n_j  - \bar t_i - \underline{{h_j\over 2} n_j r_i}=0\quad \hbox{on } \Gamma _q \quad i=1,n_b \quad j=1,n_d </math>
286
|}
287
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
288
|}
289
290
where <math display="inline">q_{ij}</math> are the generalized “fluxes” (such as the heat fluxes in a thermal  problem or the stresses in solid or fluid mechanics), <math display="inline">\bar t_i</math> are the prescribed  boundary fluxes and <math display="inline">n_j</math> are the components of the outward normal to the Neumann boundary <math display="inline">\Gamma _q</math>.
291
292
In above equations  we have underlined once more the terms introduced by the FIC approach which are essential for deriving stabilized numerical formulations.
293
294
==3 WHAT DO THE CHARACTERISTIC PARAMETERS MEAN?==
295
296
The characteristic parameters in the space and time dimension (<math display="inline">h_i</math> and <math display="inline">\delta </math>) can be interpreted as free ''intrinsic parameters'' giving the “exact” expression of the balance equations (up to first order terms) in a space-time domain of finite size.
297
298
Let us consider now the discretized solution of the modified governing equations. For simplicity we will focuss here in the simplest scalar 1D problem. The variable <math display="inline">\phi </math> is approximated as <math display="inline">\phi \simeq \hat \phi </math> where <math display="inline">\hat \phi </math> denotes the approximated solution. The values of <math display="inline">\hat \phi </math> are now expressed in terms of a finite set of parameters using any discretization procedure (finite elements, finite diferences, etc.). The discretized FIC governing equation would read now
299
300
<span id='eq-16'></span>
301
{| class="formulaSCP" style="width: 100%; text-align: left;" 
302
|-
303
| 
304
{| style="text-align: left; margin:auto;width: 100%;" 
305
|-
306
| style="text-align: center;" | <math>\hat r- \frac{h}{2}{\partial \hat r \over \partial x}-\frac{\delta }{2}{\partial \hat r \over \partial t}=r_{\bar \Omega }\quad , \qquad \hbox{in } \bar \Omega  </math>
307
|}
308
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
309
|}
310
311
where <math display="inline">\hat r:=r (\hat \phi )</math> and <math display="inline">r_{\bar \Omega }</math> is the residual of the discretized FIC equations.
312
313
The meaning of the characteristic parameters <math display="inline">h</math> and <math display="inline">\delta </math> in the discretized equation ([[#eq-16|16]]) changes  (although the same simbols than in Eqs.([[#eq-11|11]]) and ([[#eq-12|12]])  have been kept). The <math display="inline">h</math> and <math display="inline">\delta </math> parameters become now ''of the order of magnitude of the discrete domain'' where the balance laws are satisfied. In practice, this means that
314
315
<span id='eq-17'></span>
316
{| class="formulaSCP" style="width: 100%; text-align: left;" 
317
|-
318
| 
319
{| style="text-align: left; margin:auto;width: 100%;" 
320
|-
321
| style="text-align: center;" | <math>\begin{array}{c}-l \le h\le l\\ 0\le \delta \le \Delta t\end{array} </math>
322
|}
323
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
324
|}
325
326
where <math display="inline">l</math> is the grid dimension in the space discretization (i.e. the element size in the FEM or the cell size in the FDM) and <math display="inline">\Delta t</math> is the time step used to solve the transient problem.
327
328
The values of the characteristic parameters  can be found now in order to obtain an ''a correct numerical solution''. The meaning of “correct solution” must be obviously properly defined. Ideally, this would be a solution giving “exact” values at a discrete number of points (i.e. the nodes in a FE mesh). This is infortunatelly impossible for practical problems (with the exception of very simple 1D and 2D cases) and, in practice, <math display="inline">h</math> and <math display="inline">\delta </math> are computed making use of some chosen optimality rule, such as ensuring that the error of the numerical solution diminishes for appropriate values of the characteristic parameters <span id='citeF-22'></span>[[#cite-22|[22]]-<span id='citeF-27'></span>[[#cite-22|27]]]. This suffices in practice to obtain physically sound (stable) numerical results for any range of the physical parameters of the problem, always within the limitation of the discretization method chosen.
329
330
It is quite usual to accept that the characteristic parameters are ''constant'' within each element. This assumption is not justificable “a priori” and, in general, we should regard those parameters as functions of the space and time dimension and of the solution at each point of the analysis domain.
331
332
==4 FINITE ELEMENT DISCRETIZATION OF THE FIC EQUATIONS FOR ADVECTIVE-DIFFUSIVE PROBLEMS==
333
334
Let us consider the FIC governing equations for the steady-state advective-diffusive problem defined in vector form for the multidimensional case as
335
336
<span id='eq-18'></span>
337
{| class="formulaSCP" style="width: 100%; text-align: left;" 
338
|-
339
| 
340
{| style="text-align: left; margin:auto;width: 100%;" 
341
|-
342
| style="text-align: center;" | <math>r - {1\over 2}{\boldsymbol h}^T {\boldsymbol \nabla }r=0\quad \hbox{in } \Omega  </math>
343
|}
344
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
345
|}
346
347
with
348
349
<span id='eq-20a'></span>
350
{| class="formulaSCP" style="width: 100%; text-align: left;" 
351
|-
352
| 
353
{| style="text-align: left; margin:auto;width: 100%;" 
354
|-
355
| style="text-align: center;" | <math>\phi -\bar \phi =0 \quad \hbox{on } \Gamma _\phi </math>
356
|}
357
| style="width: 5px;text-align: right;white-space: nowrap;" | (20a)
358
|}
359
360
<span id='eq-20b'></span>
361
{| class="formulaSCP" style="width: 100%; text-align: left;" 
362
|-
363
| 
364
{| style="text-align: left; margin:auto;width: 100%;" 
365
|-
366
| style="text-align: center;" | <math> {\boldsymbol n}^T {\boldsymbol D} {\boldsymbol \nabla }\phi + \bar {\boldsymbol q} -  {1\over 2}{\boldsymbol h}^T{\boldsymbol n} r  =0 \quad \hbox{on } \Gamma _q</math>
367
|}
368
| style="width: 5px;text-align: right;white-space: nowrap;" | (20b)
369
|}
370
371
with
372
373
<span id='eq-21'></span>
374
{| class="formulaSCP" style="width: 100%; text-align: left;" 
375
|-
376
| 
377
{| style="text-align: left; margin:auto;width: 100%;" 
378
|-
379
| style="text-align: center;" | <math>r:= - {\boldsymbol v}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^{\!\! T} {\boldsymbol D} {\boldsymbol \nabla }\phi + Q </math>
380
|}
381
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
382
|}
383
384
In above equations <math display="inline">\boldsymbol h</math> is the characteristic length vector, <math display="inline">{\boldsymbol \nabla }</math> is the gradient vector, <math display="inline">{\boldsymbol D}</math> is the diffusivity matrix, <math display="inline">{\boldsymbol n}</math> is the normal vector and <math display="inline">{\boldsymbol v}</math> is the velocity vector.
385
386
A finite element interpolation of the unknown <math display="inline">\phi </math> can be written as
387
388
<span id='eq-22'></span>
389
{| class="formulaSCP" style="width: 100%; text-align: left;" 
390
|-
391
| 
392
{| style="text-align: left; margin:auto;width: 100%;" 
393
|-
394
| style="text-align: center;" | <math>\phi \simeq \hat \phi =\sum N_i \hat \phi _i </math>
395
|}
396
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
397
|}
398
399
where <math display="inline">N_i</math> are the shape functions and <math display="inline">\hat \phi _i</math> are the nodal values of the approximate function <math display="inline">\hat \phi </math> <span id='citeF-1'></span>[[#cite-1|[1]]].
400
401
Application of the Galerkin FE method to Equations ([[#eq-12|12]])-([[#eq-13|13]]) gives, after integrating by parts the term <math display="inline">{\boldsymbol \nabla }r</math>  (and neglecting the space derivatives of <math display="inline">{\boldsymbol h}</math>)
402
403
<span id='eq-23'></span>
404
{| class="formulaSCP" style="width: 100%; text-align: left;" 
405
|-
406
| 
407
{| style="text-align: left; margin:auto;width: 100%;" 
408
|-
409
| style="text-align: center;" | <math>\int _\Omega \!\!N_i \hat r d\Omega -\!\!\!\int _{\Gamma _q} \!\!N_i ({\boldsymbol n}^T {\boldsymbol D} {\boldsymbol \nabla }\hat \phi +\bar q_n)d\Gamma + \!\!\sum \limits _e {1\over 2}  \int _{\Omega ^e}\!\! ({\boldsymbol \nabla } {\boldsymbol h}^T  {\boldsymbol \nabla } N_i  +  N_i{\boldsymbol \nabla }^{\!\! T} {\boldsymbol h}) \hat r d\Omega =0 </math>
410
|}
411
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
412
|}
413
414
The last integral in Equation ([[#eq-23|23]]) has been expressed as the sum of the element contributions to allow for interelement discontinuities in the term <math display="inline">{\boldsymbol \nabla }\hat r</math>, where <math display="inline">\hat r = r (\hat \phi )</math> is the residual of the FE solution of the infinitesimal equations.
415
416
Note that the residual terms have disappeared from the Neumann boundary <math display="inline">\Gamma _q</math>. This is due to the consistency of the FIC terms in Equation ([[#eq-20b|20b]]).                   
417
418
The last term in the third integral involving the derivatives of <math display="inline"> \boldsymbol h</math> vanishes if <math display="inline"> \boldsymbol h</math> is assumed to be constant wihtin each element.
419
420
===4.1 Equivalence with SUPG form===
421
422
The definition of vector <math display="inline"> \boldsymbol h</math> is a crucial step as the quality of the stabilized solution depends on the module and direction of <math display="inline"> \boldsymbol h</math>.
423
424
We could, for instance,  assume that  vector <math display="inline">\boldsymbol h</math> is ''parallel'' to the velocity <math display="inline">\boldsymbol v</math>, i.e. <math display="inline">{\boldsymbol h}= h{{\boldsymbol v}\over \vert {\boldsymbol v}\vert }</math> where <math display="inline">h</math> is a characteristic length. Under these conditions, Equation ([[#eq-23|23]]) reads (for <math display="inline"> \boldsymbol h</math> assumed to be constant within each element)
425
426
<span id='eq-24'></span>
427
{| class="formulaSCP" style="width: 100%; text-align: left;" 
428
|-
429
| 
430
{| style="text-align: left; margin:auto;width: 100%;" 
431
|-
432
| style="text-align: center;" | <math>\int _\Omega N_i \hat r d\Omega - \int _{\Gamma _q} N_i ({\boldsymbol n}^T {\boldsymbol D} {\boldsymbol \nabla }\hat \phi +\bar q_n)d\Omega{+} \sum \limits _e \int _{\Omega ^e} {h\over 2\vert {\boldsymbol v}\vert  }   {\boldsymbol v}^T{\boldsymbol \nabla } N_i \hat rd\Omega =0 </math>
433
|}
434
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
435
|}
436
437
Equation ([[#eq-24|24]]) coincides precisely with the so called Streamline-Upwind-Petrov-Galerkin (SUPG) method <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-6'></span>[[#cite-6|6]],<span id='citeF-7'></span>[[#cite-7|7]],<span id='citeF-10'></span>[[#cite-10|10]]]. The ratio <math display="inline">\displaystyle{h\over 2\vert {\boldsymbol v}\vert  }</math> has dimensions of time and it is termed element ''intrinsic time'' parameter <math display="inline">\tau </math>.
438
439
It is important to note that the SUPG expression is a ''particular case'' of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction. In general, the adequate direction of <math display="inline">{\boldsymbol h}</math> is not coincident with that of <math display="inline">{\boldsymbol v}</math> and the components of <math display="inline">{\boldsymbol h}</math> introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-34'></span>[[#cite-34|34]],<span id='citeF-35'></span>[[#cite-35|35]]].
440
441
==5 INTERPRETATION OF THE DISCRETE SOLUTION OF THE FIC EQUATIONS==
442
443
Let us consider the solution of a physical problem in a space domain <math display="inline">\Omega </math>, governed by a differential equation <math display="inline">r(\phi ) =0</math> in <math display="inline">\Omega </math> with the corresponding boundary conditions. The “exact” (analytical) solution of the problem will be a function giving the sought distribution of <math display="inline">\phi </math> for any value of the geometrical and physical parameters of the problem. Obviously, since the analytical solution is difficult to find (practically impossible for real situations), an approximate numerical solution is found <math display="inline">\phi \simeq \hat \phi </math> by solving the problem <math display="inline">\hat r =0</math>, with <math display="inline">\hat r = r(\hat \phi )</math>, using a particular discretization method (such as the FEM). The distribution of <math display="inline">\phi </math> in <math display="inline">\Omega </math> is now obtained for specific values of the geometrical and physical parameters. The accuracy of the numerical solution  depends on the discretization parameters, such as the number of elements and the approximating functions chosen in the FEM. Figure [[#img-1|1]] shows a schematic representation of the distribution of <math display="inline">\hat \phi </math> along a line for different discretizations <math display="inline">M_1,M_2,\cdots , M_n</math> where <math display="inline">M_1</math> and <math display="inline">M_n</math> are the coarser and finer meshes, respectively. Obviously for <math display="inline">n</math> being sufficiently large a good approximation of <math display="inline">\phi </math> will be obtained and for <math display="inline">M_\infty </math> the numerical solution <math display="inline">\hat \phi </math> will coincide with the “exact” (and probably unreachable) analytical solution <math display="inline">\phi </math> at all points. Indeed in some problems the <math display="inline">M_\infty </math> solution can be found by a “clever” choice of the discretization parameters.
444
445
An unstable solution will occur when for some (typically coarse) discretizations, the numerical solution provides non physical or very unaccurate values of <math display="inline"> \hat \phi </math>. A situation of this kind is represented by curves <math display="inline">M_1</math> and <math display="inline">M_2</math> of the left hand side of Figure [[#img-1|1]]. These unstabilities will disappear by an appropriate mesh refinement  (curves <math display="inline">M_3,M_4 \cdots </math> in Figure [[#img-1|1]]) at the obvious increase of the computational cost.
446
447
In the FIC formulation the starting point are the modified differential equations of the problem  as previously described. These equations are however not useful to find an analytical solution, <math display="inline">\phi (x)</math>, for the physical problem. Nevertheless, the numerical solution of the FIC equation can be readily found. Moreover, by adequately choosing the values of the characteristic length parameter <math display="inline">h</math>, the numerical solution of the FIC equations will be always stable (physically sound) for any discretization level chosen.
448
449
This process is schematically represented in Figure [[#img-1|1]] where it is shown that the numerical oscillations  for the coarser discretizations <math display="inline">M_1</math> and <math display="inline">M_2</math> dissapear when using the FIC procedure.
450
451
We can  conclude the FIC approach allows us to obtain a ''better numerical solution for a given discretization''. Indeed, as in the standard infinitesimal case, the choice of <math display="inline">M_\infty </math> will yield the (unreachable) exact analytical solution and this ensures the consistency  of the method.
452
453
<div id='img-3'></div>
454
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
455
|-
456
|[[Image:Draft_Samper_674218428-Fig13con215.png|390px|Schematic representation of the numerical solution of a physical problem using standard infinitesimal calculus and finite calculus.]]
457
|- style="text-align: center; font-size: 75%;"
458
| colspan="1" | '''Figure 3:''' Schematic representation of the numerical solution of a physical problem using standard infinitesimal calculus and finite calculus.
459
|}
460
461
==6 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR==
462
463
The computation of the characteristic lengths is a crucial step as its value affect to the stability (and accuracy) of the numerical solution. This problem is common to all stabilized FE methods and different approaches to compute the stabilization parameters using typically extensions of the optimal values for simple 1D case (giving a nodally exact solution) have been proposed <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]],<span id='citeF-6'></span>[[#cite-6|6]]-<span id='citeF-21'></span>[[#cite-21|21]]].
464
465
The computation the characteristic length values in 2D and 3D problems is of much higher complexity than in 1D problems. Indeed in 1D situations <math display="inline">h</math> is an escalar (either positive or negative) while it becomes a vector 2D/3D problems. Moreover, numerical experiments indicate that in 2D/3D situations the direction of vector <math display="inline">h</math> is crucial in order to obtain stabilized numerical solutions in problems where boundary layers and arbitrary internal sharp layers  exist. It is well known that the SUPG assumption (<math display="inline">h</math> being parallel to <math display="inline">u</math>) generally does not suffice to give stable results in those cases <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]],<span id='citeF-8'></span>[[#cite-8|8]],<span id='citeF-10'></span>[[#cite-10|10]],<span id='citeF-34'></span>[[#cite-34|34]],<span id='citeF-35'></span>[[#cite-35|35]]]. Conversely, the numerical results in those cases are more insensitive to the module of <math display="inline">h</math> which can be taken of the order of a typical element dimension.
466
467
Much effort has been spent by the author in the past in order to derive numerical schemes for computing vector <math display="inline">h</math> in an iterative manner. The interested reader can find details on the different procedures in <span id='citeF-22'></span>[[#cite-22|[22]]-<span id='citeF-30'></span>[[#cite-30|30]]].
468
469
We present here a new approach for computing vector <math display="inline">h</math> which is general and applicable to all FIC equations in mechanics.
470
471
The basis of the method is to ''assume that vector <math display="inline">h</math> is a function of the gradient of the numerical solution''. The simplest choice for  convection-diffusion problems is
472
<span id='eq-25'></span>
473
{| class="formulaSCP" style="width: 100%; text-align: left;" 
474
|-
475
| 
476
{| style="text-align: left; margin:auto;width: 100%;" 
477
|-
478
| style="text-align: center;" | <math>{\boldsymbol h}={\boldsymbol H} {\boldsymbol \nabla } \hat \phi  </math>
479
|}
480
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
481
|}
482
483
where <math display="inline">{\boldsymbol H}</math> is a matrix which terms are a function of the element size and the inverse of the gradient vector components. The following simple expression of <math display="inline">{\boldsymbol H}</math> for 2D problems has found to give excellent numerical results in practice <span id='citeF-43'></span>[[#cite-43|[43]]]
484
<span id='eq-26'></span>
485
{| class="formulaSCP" style="width: 100%; text-align: left;" 
486
|-
487
| 
488
{| style="text-align: left; margin:auto;width: 100%;" 
489
|-
490
| style="text-align: center;" | <math>{\boldsymbol H} =h\left[\begin{matrix}\displaystyle \left({\partial \hat \phi \over \partial x}\right)^{-1} &0\\ 0 & \displaystyle \left({\partial \hat \phi \over \partial y}\right)^{-1}\\\end{matrix}\right] </math>
491
|}
492
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
493
|}
494
495
where <math display="inline">h</math> is a characteristic length parameter. As mentioned above the value of this length is not so relevant in practice. The rationale for choosing Eqs.([[#eq-25|25]]) is explained in <span id='citeF-43'></span>[[#cite-43|[43]]].
496
497
A simple expression for computing the length <math display="inline">h^e</math> for each element is
498
<span id='eq-27'></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>h^e=\max \left|{\boldsymbol l}_i^T {{\boldsymbol \nabla } \hat \phi \over \vert {\boldsymbol \nabla } \hat \phi \vert }\right|\quad ,\qquad i=1,n_l </math>
505
|}
506
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
507
|}
508
509
where <math display="inline">{\boldsymbol l}_i</math> are the element side vectors and <math display="inline">n_l</math> is the number of sides for each element (i.e. <math display="inline">n_l=3</math> for triangles, etc.).
510
511
Substituting Eqs.([[#eq-25|25]]) into the weighted residual form Eqs.([[#eq-23|23]]) gives
512
<span id='eq-28'></span>
513
{| class="formulaSCP" style="width: 100%; text-align: left;" 
514
|-
515
| 
516
{| style="text-align: left; margin:auto;width: 100%;" 
517
|-
518
| style="text-align: center;" | <math>\int _\Omega N_i \hat r d\Omega + \sum \limits _e {1\over 2} \int _{\Omega ^e} h \left[{\boldsymbol \nabla }^{\!\! T} N_i {\boldsymbol H} {\boldsymbol \nabla }\hat \phi + N_i {\boldsymbol \nabla }^{\!\! T} {\boldsymbol H} {\boldsymbol \nabla }\hat \phi \right]\hat r d\Omega + b.t.=0 </math>
519
|}
520
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
521
|}
522
523
where ''b.t.'' stands for the boundary  integral terms.
524
525
Let us assume now that a linear interpolation is taken for <math display="inline">\phi </math> in Eq.([[#eq-22|22]]). This allows to neglect the second term of the second integral in ([[#eq-28|28]]). The new expression can be written as
526
<span id='eq-29'></span>
527
{| class="formulaSCP" style="width: 100%; text-align: left;" 
528
|-
529
| 
530
{| style="text-align: left; margin:auto;width: 100%;" 
531
|-
532
| style="text-align: center;" | <math>\int _\Omega N_i \hat r d\Omega + \sum \limits _e {1\over 2} \int _{\Omega ^e} ({\boldsymbol \nabla }^{\!\! T} N_i)  \hat r {\boldsymbol H}   {\boldsymbol \nabla }\phi d\Omega + b.t.=0 </math>
533
|}
534
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
535
|}
536
537
We see now clearly that the second integral has the form of a Laplacian matrix where the term <math display="inline">\displaystyle{\hat r  \over 2}{\boldsymbol H}</math> takes the role of a diffusivity matrix.
538
539
Let us integrate now the diffusion term within the first integral of Eq.([[#eq-29|29]]). This gives after small algebra
540
<span id='eq-30'></span>
541
{| class="formulaSCP" style="width: 100%; text-align: left;" 
542
|-
543
| 
544
{| style="text-align: left; margin:auto;width: 100%;" 
545
|-
546
| style="text-align: center;" | <math>\sum \limits _e \int _{\Omega ^e} \left[N_i {\boldsymbol v}^T {\boldsymbol \nabla }\hat \phi + ({\boldsymbol \nabla }^{\!\! T} N_i) {\boldsymbol D}^* {\boldsymbol \nabla } \hat \phi \right]d\Omega - \int _\Omega N_i Q d\Omega + \int _{\Gamma _q} N_i \bar q_n d\Gamma =0 </math>
547
|}
548
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
549
|}
550
551
where the diffusivity matrix <math display="inline">{\boldsymbol D}^*</math> is
552
<span id='eq-31'></span>
553
{| class="formulaSCP" style="width: 100%; text-align: left;" 
554
|-
555
| 
556
{| style="text-align: left; margin:auto;width: 100%;" 
557
|-
558
| style="text-align: center;" | <math>{\boldsymbol D}^* = {\boldsymbol D} + {\vert \hat r \vert \over 2}{\boldsymbol H}  </math>
559
|}
560
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
561
|}
562
563
and <math>I</math> is the unit matrix. In Eq.([[#eq-29|29]]) the absolute value of <math display="inline">\vert  \hat r \vert </math> is taken to ensure a positive value of the new diffusivity terms.
564
565
Eq.([[#eq-30|30]]) yields the final system of discretized equations in the standard form
566
<span id='eq-32'></span>
567
{| class="formulaSCP" style="width: 100%; text-align: left;" 
568
|-
569
| 
570
{| style="text-align: left; margin:auto;width: 100%;" 
571
|-
572
| style="text-align: center;" | <math>{\boldsymbol K} {\boldsymbol a} ={\boldsymbol f} </math>
573
|}
574
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
575
|}
576
577
where the stiffness matrix <math>K</math> and the nodal vector <math>f</math> are assembled from the element contributions
578
<span id='eq-33'></span><span id='eq-34'></span>
579
{| class="formulaSCP" style="width: 100%; text-align: left;" 
580
|-
581
| 
582
{| style="text-align: left; margin:auto;width: 100%;" 
583
|-
584
| style="text-align: center;" | <math>{\boldsymbol K}_{ij}^e = \int _{\Omega ^e} ({\boldsymbol \nabla }^{\!\! T} N_i) {\boldsymbol D}^* {\boldsymbol \nabla }N_jd\Omega </math>
585
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
586
|-
587
| style="text-align: center;" | <math> {\boldsymbol f}_{i}^e = \int _{\Omega ^e}N_i Q d\Omega  </math>
588
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
589
|}
590
|}
591
592
An iterative algorithm giving a stabilized solution can be implemented as follows
593
594
<ol>
595
596
<li>Compute an initial solution <math display="inline">{\boldsymbol a}^1</math> solving Eq.([[#eq-33|33]]) for an initial value of <math display="inline">{\boldsymbol D}^*</math>. Typically one can chose
597
<span id='eq-35'></span>
598
{| class="formulaSCP" style="width: 100%; text-align: left;" 
599
|-
600
| 
601
{| style="text-align: left; margin:auto;width: 100%;" 
602
|-
603
| style="text-align: center;" | <math>
604
605
^1{\boldsymbol D}^*={\boldsymbol D}+{h^e \over 2}\left[\begin{matrix}\vert u\vert &0\\ 0 & \vert v\vert \\\end{matrix}\right] </math>
606
|}
607
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
608
|}</li>
609
<li>Compute enhanced values of the gradient <math display="inline">\overline{{\boldsymbol \nabla }\hat \phi }</math>. This can be performed using a derivative recovery scheme <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-38'></span>[[#cite-38|38]]]. </li>
610
611
<li>Compute updated element values of the characteristic length vector <math display="inline">{}^{i+1}{\boldsymbol h}</math>
612
<span id='eq-36'></span>
613
{| class="formulaSCP" style="width: 100%; text-align: left;" 
614
|-
615
| 
616
{| style="text-align: left; margin:auto;width: 100%;" 
617
|-
618
| style="text-align: center;" | <math>
619
620
{}^{i+1}{\boldsymbol h}^e = h^e \left.\begin{matrix}{}^{^i}\\\\\end{matrix}\right.\!\!\!\left[{\overline{{\boldsymbol \nabla }\hat \phi }\over \vert \overline{{\boldsymbol \nabla }\hat \phi }\vert }\right]^e </math>
621
|}
622
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
623
|}</li>
624
625
where <math display="inline">{}^i\bar{(\cdot )}^e</math> denotes mean enhanced values for element <math display="inline">e</math> for the <math display="inline">i</math>th iteration.
626
627
<li>Compute the element residuals using the enhanced derivative field by
628
<span id='eq-37'></span>
629
{| class="formulaSCP" style="width: 100%; text-align: left;" 
630
|-
631
| 
632
{| style="text-align: left; margin:auto;width: 100%;" 
633
|-
634
| style="text-align: center;" | <math>
635
636
{}^i \bar{\hat r}^e = \int _{\Omega ^e} \left[{}^i \bar{\hat r} - {1\over 2} [{}^{i+1}{\boldsymbol h}^e]^T {\boldsymbol \nabla } {}^i \bar{\hat r}\right]d\Omega  </math>
637
|}
638
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
639
|}</li>
640
641
<li>Check for convergence of element residuals
642
<span id='eq-38'></span>
643
{| class="formulaSCP" style="width: 100%; text-align: left;" 
644
|-
645
| 
646
{| style="text-align: left; margin:auto;width: 100%;" 
647
|-
648
| style="text-align: center;" | <math>
649
650
{\left[\sum \limits _e ({}^i \bar{\hat r}^e)^2\right]^{1/2}\over NQ \Omega ^e }\le \varepsilon _r </math>
651
|}
652
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
653
|}</li>
654
655
where <math display="inline">\varepsilon _r</math> is a prescribed tolerance (typically <math display="inline">\varepsilon _r \simeq 10^{-4}</math>) and <math display="inline">N</math> is the total number of elements in the mesh.
656
657
The term in the denominator in Eq.([[#eq-38|38]]) is chosen so as to scale the residual error. For <math display="inline">Q=0</math> the denominator should be replaced by <math display="inline">N[\Omega ^e]^{1/2} \vert {\boldsymbol v}_{max}\vert \bar \phi _{max}</math> where <math display="inline">{\boldsymbol v}_{max}</math> is the maximum value of the velocity vector in the mesh, <math display="inline">\bar \phi _{max}</math> is the maximum prescribed value of <math display="inline">\phi </math> at the Dirichlet boundary and <math display="inline">n_d=2/3</math> for 2D/3D problems.
658
659
<li>Repeat steps 1&#8211;5 until convergence is satisfied using the updated value of <math display="inline">{\boldsymbol D}^*</math>
660
<span id='eq-39'></span>
661
{| class="formulaSCP" style="width: 100%; text-align: left;" 
662
|-
663
| 
664
{| style="text-align: left; margin:auto;width: 100%;" 
665
|-
666
| style="text-align: center;" | <math>
667
668
^i{\boldsymbol D}^* = {\boldsymbol D} + {\vert {}^i \bar{\hat r} \vert \over 2} {}^i {\boldsymbol H} </math>
669
|}
670
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
671
|}</li>
672
673
</ol>
674
675
Numerical experimentals have shown that above process yields a converged stabilized solution in 2&#8211;3 iterations. Details and extensions  of this scheme including numerical results can be found in  <span id='citeF-43'></span>[[#cite-43|[43]]].
676
677
==7 GENERALIZATION OF THE FIC STABILIZATION PROCESS WITH GRADIENT ORIENTED CHARACTERISTIC LENGTH VECTORS==
678
679
The FIC stabilization process described in previous section can be generalized by choosing vector '''h''' in the direction of the solution gradient as follows.
680
681
Let us consider a more general expression of the FIC equation ([[#eq-11|11]]) for a multidimensional problem where a different expression for <math>h</math> is chosen for each balance equation (for simplicity we restrict onselves to steady state problems only)
682
<span id='eq-40'></span>
683
{| class="formulaSCP" style="width: 100%; text-align: left;" 
684
|-
685
| 
686
{| style="text-align: left; margin:auto;width: 100%;" 
687
|-
688
| style="text-align: center;" | <math>r_i - {1\over 2} {\boldsymbol \nabla }^T {\boldsymbol h}_i r_i =0\quad \hbox{no sum in }i\quad ;\quad i=1,n_b~;~j=1,n_d </math>
689
|}
690
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
691
|}
692
693
where <math display="inline">n_b</math> and <math display="inline">n_d</math> were defined in Eq.([[#eq-11|11]]).
694
695
We choose the characteristic length vectors <math display="inline">{\boldsymbol h}_i</math> as
696
<span id='eq-41'></span>
697
{| class="formulaSCP" style="width: 100%; text-align: left;" 
698
|-
699
| 
700
{| style="text-align: left; margin:auto;width: 100%;" 
701
|-
702
| style="text-align: center;" | <math>{\boldsymbol h}_i=h_i {{\boldsymbol \nabla } u_i\over \vert {\boldsymbol \nabla } u_i\vert }\qquad \hbox{no sum in }i </math>
703
|}
704
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
705
|}
706
707
where <math display="inline">u_j</math> is an appropriate variable corresponding to the <math display="inline">i</math>th balance equation. For instance in a fluid mechanics problem <math display="inline">u_i</math> would be the velocity along the <math display="inline">i</math>th axis.
708
709
The weak form of Eq.([[#eq-40|40]]) is obtained after finite element discretization as
710
<span id='eq-42'></span>
711
{| class="formulaSCP" style="width: 100%; text-align: left;" 
712
|-
713
| 
714
{| style="text-align: left; margin:auto;width: 100%;" 
715
|-
716
| style="text-align: center;" | <math>\int _{\Omega ^e}  N_k \hat r_i d\Omega + \sum \limits _e {1\over 2} \int _{\Omega ^e} {h_i\hat r_i  \over \vert {\boldsymbol \nabla }\hat u_i\vert } {\boldsymbol \nabla }^{\!\! T} N_k {\boldsymbol \nabla }\hat u_i d\Omega + b.t.=0\quad \hbox{no sum in }i </math>
717
|}
718
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
719
|}
720
721
We note that the stabilization term plays the role of a Laplacian which introduces a residual-type diffusion into the Galerkin equation.
722
723
The procedure is particularized next for the equations of an incompressible fluid.
724
725
==8 FIC METHOD FOR INCOMPRESSIBLE FLUID MECHANICS==
726
727
The FIC method can be applied to derive the modified equations of momentum,  mass and energy conservation in fluid mechanics. The general form of these equations for a compressible fluid was presented in <span id='citeF-22'></span>[[#cite-22|[22]],<span id='citeF-29'></span>[[#cite-29|29]],<span id='citeF-32'></span>[[#cite-32|32]]]. We will consider here the particular case of a ''viscous incompressible fluid''. The FIC equations for the momentum and mass balance in this case can be written as (neglecting time stabilization terms)
728
729
'''Momentum'''
730
<span id='eq-43'></span>
731
{| class="formulaSCP" style="width: 100%; text-align: left;" 
732
|-
733
| 
734
{| style="text-align: left; margin:auto;width: 100%;" 
735
|-
736
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_{m_j}^i{\partial r_{m_i} \over \partial x_j}}=0 \qquad \hbox{no sum in }i </math>
737
|}
738
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
739
|}
740
741
'''Mass balance'''
742
<span id='eq-44'></span>
743
{| class="formulaSCP" style="width: 100%; text-align: left;" 
744
|-
745
| 
746
{| style="text-align: left; margin:auto;width: 100%;" 
747
|-
748
| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_{d_j} {\partial r_d \over \partial x_j}}=0  </math>
749
|}
750
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
751
|}
752
753
where
754
<span id='eq-45'></span><span id='eq-46'></span>
755
{| class="formulaSCP" style="width: 100%; text-align: left;" 
756
|-
757
| 
758
{| style="text-align: left; margin:auto;width: 100%;" 
759
|-
760
| style="text-align: center;" | <math>r_{m_i} = \rho \left({\partial v_i \over \partial t}+v_j {\partial v_i \over \partial x_j}\right)+ {\partial p \over \partial x_i}- {\partial s_{ij} \over \partial x_j}-b_i</math>
761
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
762
|-
763
| style="text-align: center;" | <math> r_d = {\partial v_i \over \partial x_i}\qquad i,j = 1, n_d </math>
764
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
765
|}
766
|}
767
768
Above <math display="inline">v_i</math> is the velocity along the ith global axis, <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">p</math> is the absolute pressure (defined positive in compression), <math display="inline">b_i</math> are body forces and <math display="inline">s_{ij}</math> are the viscous deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression
769
<span id='eq-47'></span>
770
{| class="formulaSCP" style="width: 100%; text-align: left;" 
771
|-
772
| 
773
{| style="text-align: left; margin:auto;width: 100%;" 
774
|-
775
| style="text-align: center;" | <math>s_{ij}=2\mu \left(\dot \varepsilon _{ij} - \delta _{ij} {1\over 3} {\partial v_k \over \partial x_k}\right) </math>
776
|}
777
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
778
|}
779
780
where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are
781
<span id='eq-48'></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>\dot \varepsilon _{ij}={1\over 2} \left({\partial v_i \over \partial x_j}+{\partial v_j \over \partial x_i}\right) </math>
788
|}
789
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
790
|}
791
792
The FIC boundary conditions are written as
793
<span id='eq-49'></span>
794
{| class="formulaSCP" style="width: 100%; text-align: left;" 
795
|-
796
| 
797
{| style="text-align: left; margin:auto;width: 100%;" 
798
|-
799
| style="text-align: center;" | <math>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_{m_j} n_j r_{m_i}}=0 \quad \hbox{on }\Gamma _t </math>
800
|}
801
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
802
|}
803
<span id='eq-50'></span>
804
{| class="formulaSCP" style="width: 100%; text-align: left;" 
805
|-
806
| 
807
{| style="text-align: left; margin:auto;width: 100%;" 
808
|-
809
| style="text-align: center;" | <math>v_j - \bar v_j =0 \quad \hbox{on }\Gamma _u </math>
810
|}
811
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
812
|}
813
814
and the initial condition <math display="inline">v_j =v_j^0</math> for <math display="inline">t=t_0</math>.
815
816
In Equation ([[#eq-45|45]]) <math display="inline">\sigma _{ij}= s_{ij}-p\delta _{ij}</math> are the total stresses, <math display="inline">t_i</math> and <math display="inline">\bar u_j</math> are prescribed tractions and displacements on the boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _u</math>, respectively and <math display="inline">n_j</math> are the components of the unit normal vector to the boundary.  The sign in front the stabilization term in Equation ([[#eq-49|49]]) is positive due to the definition of <math display="inline">r_{m_i}</math> in Eq.([[#eq-45|45]]).
817
818
The <math display="inline">h_{m_j}^i</math> and <math display="inline">h_{d_j}</math> in above equations are characteristic lengths of the domain where balance of momentum and mass is enforced. In Equation ([[#eq-49|49]]) the lengths <math display="inline">h_{m_j}</math> define the domain where equilibrium of boundary tractions is established <span id='citeF-22'></span>[[#cite-22|[22]]].
819
820
Equations ([[#eq-43|43]])&#8211;([[#eq-50|50]]) are the starting  point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations using an equal order interpolation for the velocity and the pressure variables.
821
822
The weighted residual form of the momentum and mass balance equations can be written as
823
<span id='eq-51'></span> 
824
{| class="formulaSCP" style="width: 100%; text-align: left;" 
825
|-
826
| 
827
{| style="text-align: left; margin:auto;width: 100%;" 
828
|-
829
| style="text-align: center;" | <math>\int _\Omega \delta v_i \left[r_{m_i} - {h^i_{m_j}\over 2} {\partial r_{m_i} \over \partial x_j}\right]d\Omega + \int _{\Gamma _t} \delta v_i \left(\sigma _{ij}n_j - t_i + {h_j\over 2} n_j r_{m_i}\right)d\Gamma =0 </math>
830
|}
831
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
832
|}
833
<span id='eq-52'></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>\int _\Omega q \left(r_d - {h_{d_j}\over 2} {\partial r_{d} \over \partial x_j}\right)d\Omega =0 </math>
840
|}
841
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
842
|}
843
844
where <math display="inline">\delta v_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocity and virtual pressure fields. Here and in the following no sum applies to the <math display="inline">i</math>th superindex <math display="inline">i</math> of <math display="inline">h^i_{m_j}</math>.
845
846
Integrating by parts Eqs.([[#eq-51|51]]) and ([[#eq-52|52]]) gives
847
<span id='eq-53'></span>
848
{| class="formulaSCP" style="width: 100%; text-align: left;" 
849
|-
850
| 
851
{| style="text-align: left; margin:auto;width: 100%;" 
852
|-
853
| style="text-align: center;" | <math>\int _\Omega \delta v_i r_{m_i}d\Omega + \int _{\Gamma _t} \delta v_i \left(\sigma _{ij}n_j - t_i\right)d\Gamma + \sum \limits _e \int _{\Omega ^e} { h^i_{m_j}\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i}d\Omega =0 </math>
854
|}
855
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
856
|}
857
<span id='eq-54'></span>
858
{| class="formulaSCP" style="width: 100%; text-align: left;" 
859
|-
860
| 
861
{| style="text-align: left; margin:auto;width: 100%;" 
862
|-
863
| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \sum \limits _e \int _{\Omega ^e} {h_{d_j}\over 2} {\partial q \over \partial x_j}d\Omega =0 </math>
864
|}
865
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
866
|}
867
868
In the derivation of Eqs.([[#eq-53|53]]) and ([[#eq-54|54]]) the following assumptions have been made.
869
870
<ol>
871
872
<li><math display="inline">h_{m_j}^i \equiv h_{m_j}</math> at the Neumann boundary. This allows to elliminate the residual of the momentum equations at that boundary after the integration by parts. </li>
873
<li>The volumetric strain rate <math display="inline">r_d</math> vanishes at the boundary. </li>
874
<li>The characteristic lengths <math display="inline">h_{m_j}^i</math> and <math display="inline">h_{d_j}</math> are constant within each element. </li>
875
876
</ol>
877
878
===8.1 Gradient-form of the characteristic lengths===
879
880
The following gradient-based expressions are taken for the characteristic lengths in the momentum and mass balance equations
881
<span id='eq-55a'></span> <span id='eq-55b'></span>
882
{| class="formulaSCP" style="width: 100%; text-align: left;" 
883
|-
884
| 
885
{| style="text-align: left; margin:auto;width: 100%;" 
886
|-
887
| style="text-align: center;" | <math>{\boldsymbol h}^i_{m} = h^i_m {{\boldsymbol \nabla }v_i\over \vert {\boldsymbol \nabla }v_i\vert }\quad \hbox{no sum in }i</math>
888
|}
889
| style="width: 5px;text-align: right;white-space: nowrap;" | (55a)
890
|}
891
892
{| class="formulaSCP" style="width: 100%; text-align: left;" 
893
|-
894
| 
895
{| style="text-align: left; margin:auto;width: 100%;" 
896
|-
897
| style="text-align: center;" | <math>{\boldsymbol h}_{d} = h_d  {{\boldsymbol \nabla }p\over \vert {\boldsymbol \nabla }p\vert }</math>
898
|}
899
| style="width: 5px;text-align: right;white-space: nowrap;" | (55b)
900
|}
901
902
The distances <math display="inline">h^i_{m}</math> and <math display="inline">h_d</math> are computed as
903
<span id='eq-56a'></span> <span id='eq-56b'></span>
904
{| class="formulaSCP" style="width: 100%; text-align: left;" 
905
|-
906
| 
907
{| style="text-align: left; margin:auto;width: 100%;" 
908
|-
909
| style="text-align: center;" | <math>h^i_{m}= \hbox{max} \left|{\boldsymbol l}_j^T {{\boldsymbol \nabla } v_i\over \vert{\boldsymbol \nabla } v_i\vert }  \right|\quad j=1,n_l</math>
910
|}
911
| style="width: 5px;text-align: right;white-space: nowrap;" | (56a)
912
|}
913
914
{| class="formulaSCP" style="width: 100%; text-align: left;" 
915
|-
916
| 
917
{| style="text-align: left; margin:auto;width: 100%;" 
918
|-
919
| style="text-align: center;" | <math>h_d = \hbox{max} \left|{\boldsymbol l}_j^T {{\boldsymbol \nabla } p\over \vert {\boldsymbol \nabla } p\vert } \right|\quad j=1,n_l </math>
920
|}
921
| style="width: 5px;text-align: right;white-space: nowrap;" | (56b)
922
|}
923
924
Substituting Eqs.([[#eq-55|55]]) into ([[#eq-53|53]]) and ([[#eq-54|54]]) gives after integration by parts of the deviatoric stress and pressure term of <math display="inline">r_{m_i}</math> in Eq.([[#eq-53|53]])
925
<span id='eq-57a'></span> <span id='eq-57b'></span>
926
{| class="formulaSCP" style="width: 100%; text-align: left;" 
927
|-
928
| 
929
{| style="text-align: left; margin:auto;width: 100%;" 
930
|-
931
| style="text-align: center;" | <math> \begin{array}{r} \displaystyle \int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ \delta \dot \varepsilon _{ij}(\tau _{ij}- \delta _{ij}p )\right]\, d\Omega -  \int _{\Omega } \delta v_i b_i \,d\Omega -\\ \displaystyle - \int _{\Gamma _t} \delta v_i t_i\,d\Gamma + \sum \limits _e \!\int _{\Omega ^e} {h_m^i r_{m_i}\over 2\vert {\boldsymbol \nabla } v_i\vert }{\partial \delta v_i \over \partial x_j} {\partial v_i \over \partial x_j} \,d\Omega =0 \end{array} </math>
932
|}
933
| style="width: 5px;text-align: right;white-space: nowrap;" | (57a)
934
|}
935
936
{| class="formulaSCP" style="width: 100%; text-align: left;" 
937
|-
938
| 
939
{| style="text-align: left; margin:auto;width: 100%;" 
940
|-
941
| style="text-align: center;" | <math>\int _\Omega q r_d \,d\Omega + \sum \limits _e \int _{\Omega ^e} {h_d r_d\over 2\vert {\boldsymbol \nabla } p\vert }{\partial q \over \partial x_j} {\partial p \over \partial x_j} \,d\Omega =0</math>
942
|}
943
| style="width: 5px;text-align: right;white-space: nowrap;" | (57b)
944
|}
945
946
We see clearly that the stabilization terms in both equations take the form of a Laplacian matrix as it is desirable.
947
948
We introduce now a standard equal order linear finite element interpolation of the velocity and pressure fields as
949
<span id='eq-58'></span>
950
{| class="formulaSCP" style="width: 100%; text-align: left;" 
951
|-
952
| 
953
{| style="text-align: left; margin:auto;width: 100%;" 
954
|-
955
| style="text-align: center;" | <math>v_i = \sum \limits _{j=1}^{n} N_j \bar v_i^j\quad ,\quad p = \sum \limits _{j=1}^{n} N_j \bar p_j </math>
956
|}
957
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
958
|}
959
960
where <math display="inline">N_j</math> are the linear shape functions, <math display="inline">n</math> is the number of nodes per element and <math display="inline">\overline{(\cdot )}</math> denotes nodal variables.
961
962
Substituting Eqs.([[#eq-58|58]]) into ([[#eq-57|57]]) leads to the following system of equations
963
<span id='eq-59a'></span> <span id='eq-59b'></span>
964
{| class="formulaSCP" style="width: 100%; text-align: left;" 
965
|-
966
| 
967
{| style="text-align: left; margin:auto;width: 100%;" 
968
|-
969
| style="text-align: center;" | <math>{\boldsymbol M} \dot{\bar{\boldsymbol v}} + [{\boldsymbol K}(\mu ) + {\boldsymbol A} (\bar{\boldsymbol v})+\bar {\boldsymbol L} (\bar {\boldsymbol v} )] \bar{\boldsymbol v}- {\boldsymbol G} \bar {\boldsymbol p} ={\boldsymbol f} </math>
970
|}
971
| style="width: 5px;text-align: right;white-space: nowrap;" | (59a)
972
|}
973
974
{| class="formulaSCP" style="width: 100%; text-align: left;" 
975
|-
976
| 
977
{| style="text-align: left; margin:auto;width: 100%;" 
978
|-
979
| style="text-align: center;" | <math>{\boldsymbol G}^T \bar{\boldsymbol v}+ {\boldsymbol L} (\tau _p ) \bar {\boldsymbol p} = 0</math>
980
|}
981
| style="width: 5px;text-align: right;white-space: nowrap;" | (59b)
982
|}
983
984
where for 2D problems
985
<span id='eq-60'></span>
986
{| class="formulaSCP" style="width: 100%; text-align: left;" 
987
|-
988
| 
989
{| style="text-align: left; margin:auto;width: 100%;" 
990
|-
991
| style="text-align: center;" | <math>\displaystyle M_{ij}\!\!\!\!=\!\!\!\!  \int _{\Omega ^e} \rho N_iN_j d\Omega \quad ,\quad A_{ij}= \int _{\Omega ^e} N_i \rho {\boldsymbol v}^T {\boldsymbol \nabla } N_j d\Omega \quad ,\quad  {\boldsymbol \nabla }= \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2}\right]^T </math>
992
|-
993
| style="text-align: center;" | 
994
|-
995
| style="text-align: center;" | <math> \bar{\boldsymbol L}_{ij}\!\!\!\!=\!\!\!\! \left[\begin{matrix}{\boldsymbol L}_{ij}^1 &0&0\\ 0&{\boldsymbol L}_{ij}^2 &0\\ 0&0& {\boldsymbol L}_{ij}^3\\\end{matrix}\right] d\Omega \quad ,\quad L_{ij}^k=\int _{\Omega ^e} \tau _k ({\boldsymbol \nabla }^{\!\! T} N_i) {\boldsymbol \nabla }N_jd\Omega \quad ,\quad \tau _i={\vert h_m^i r_{m_i}\vert \over 2\vert {\boldsymbol \nabla } v_i\vert }</math>
996
|-
997
| style="text-align: center;" | 
998
|-
999
| style="text-align: center;" | <math> \displaystyle {K}_{ij}\!\!\!\!=\!\!\!\! \int _{\Omega ^e}\mu ({\boldsymbol \nabla }^{\!\! T}N_i) {\boldsymbol \nabla } N_j d\Omega  \quad ,\quad \displaystyle {\boldsymbol G}_{ij}= \int _{\Omega ^e}({\boldsymbol \nabla } N_i)N_j d\Omega  </math>
1000
|-
1001
| style="text-align: center;" | 
1002
|-
1003
| style="text-align: center;" | <math> \displaystyle L_{ij}\!\!\!\!=\!\!\!\!\int _{\Omega ^e} \tau _p ({\boldsymbol \nabla }^{\!\! T} N_i) {\boldsymbol \nabla } N_j d\Omega \quad ,\quad \tau _p= {\vert h_d r_d\vert \over 2\vert{\boldsymbol \nabla }p\vert } </math>
1004
|-
1005
| style="text-align: center;" | 
1006
|-
1007
| style="text-align: center;" | <math> \displaystyle {\boldsymbol f}_i\!\!\!\! =\!\!\!\! \int _{\Omega ^e} N_i {\boldsymbol b}d\Omega + \int _{\Gamma ^e}N_i {\boldsymbol t} d\Gamma \quad ,\quad  {\boldsymbol b}=[b_1,b_2]^T \quad ,\quad {\boldsymbol t}=[t_1,t_2]^T  </math>
1008
|}
1009
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
1010
|}
1011
1012
It is interesting to analyze the steady-state form of Eqs.([[#eq-59|59]]) for the Stokes flow case where the convective terms are neglected. Now  the convective matrix <math>A</math> and the stabilization matrix <math>L</math> can be made equal to zero in the momentum equations and the resulting system can be written as
1013
<span id='eq-61'></span>
1014
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1015
|-
1016
| 
1017
{| style="text-align: left; margin:auto;width: 100%;" 
1018
|-
1019
| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K} (\mu ) & -{\boldsymbol G}\\ -{\boldsymbol G}^T & - {\boldsymbol L} (\tau _p)\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol h}\\ \bar {\boldsymbol p}\\\end{matrix}\right\}= \left\{\begin{matrix}\bar {\boldsymbol f}\\ {\boldsymbol 0}\\\end{matrix}\right\} </math>
1020
|}
1021
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
1022
|}
1023
1024
The stability of the numerical solution is ensured by the presence of matrix <math>L</math> guaranteeing a positive definitiveness of the equation system ''for any choice of the approximations for v and'' <math display="inline">p</math>, thus overcoming the Babuska-Brezzi conditions <span id='citeF-1'></span>[[#cite-1|[1]]].
1025
1026
The extension of above procedure to derive stabilized equations for incompressible solid mechanics problems is straight forward making use of the well known analogy between the Stokes equations for an incompressible flow and those of incompressible elasticity.
1027
1028
Applications of the FIC method here proposed to incompressible problems in fluid and solid mechanics can be found in <span id='citeF-44'></span>[[#cite-44|[44]],<span id='citeF-45'></span>[[#cite-45|45]]].
1029
1030
==9 CONCLUDING REMARKS==
1031
1032
We have presented in this paper the possibilities of the finite calculus (FIC) method for deriving stabilized finite element formulations for a variety of problems in mechanics. In all cases the modified differential equations derived via the FIC method yield naturally a discretized system of equations with intrinsic stabilizations properties. These equations are more general than other standard stabilization methods (such as SUPG) and they allow to obtain correct numerical solutions for complex problems involving boundary layers and sharp internal layers. The key to the succes of the FIC method is the correct selection of the characteristic vector '''h'''. We have presented a new gradient-based definition of '''h''' which introduces naturally stabilization terms of laplacian-type in the equations for convective-diffusive transport and incompressible fluid flow problems.
1033
1034
Numerical examples with applications of the formulation here presented can be found in <span id='citeF-43'></span>[[#cite-43|43]].
1035
1036
==ACKNOWLEDGEMENTS==
1037
1038
Thanks are given to Profs. J. García, S.R. Idelsohn, R.L. Taylor and O.C. Zienkiewicz for many useful discussions.
1039
1040
==References==
1041
1042
<div id="cite-1"></div>
1043
[[#citeF-1|[1]]] Zienkiewicz OC, Taylor RL. ''The finite element method''. 5th Edition, 3 Volumes, Butterworth&#8211;Heinemann, (2000).
1044
1045
<div id="cite-2"></div>
1046
[[#citeF-2|[2]]]  Hirsch C. ''Numerical computation of internal and external flow'', J. Wiley, Vol. 1 1988, Vol. 2, 1990.
1047
1048
<div id="cite-3"></div>
1049
[[#citeF-3|[3]]]  Oñate E, Idelsohn SR, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Applications to convective transport and fluid flow. ''Int. J. Num. Meth. Engng.'', 39, 3839&#8211;3866, (1996).
1050
1051
<div id="cite-4"></div>
1052
[[#citeF-4|[4]]]  Oñate E, Idelsohn S. A mesh free finite point method for advective-diffusive transport and fluid flow problems. ''Comput. Mechanics'', 21, 283&#8211;292, (1988).
1053
1054
<div id="cite-5"></div>
1055
[[#citeF-5|[5]]]  Oñate E, Sacco C, Idelsohn S. A finite point method for incompressible flow problems. ''Computing and Visualization in Science'', 2, 67&#8211;75, (2000).
1056
1057
<div id="cite-6"></div>
1058
[[#citeF-6|[6]]] Brooks A, Hughes TJR. Streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis  on the incompressible Navier-Stokes equations. ''Comput. Methods Appl.  Mech. Engrg.'', 32, 199&#8211;259, (1982).
1059
1060
<div id="cite-7"></div>
1061
[[#citeF-7|[7]]]  Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', 156, 185&#8211;210, (1998).
1062
1063
<div id="cite-8"></div>
1064
[[#citeF-8|[8]]] Hughes TJR, Mallet M. A new finite element formulations  for computational fluid dynamics: III. The generalized streamline operator  for multidimensional advective-diffusive systems. ''Comput Methods Appl.  Mech. Engrg.'', 58, 305&#8211;328, (1986a).
1065
1066
<div id="cite-9"></div>
1067
[[#citeF-9|[9]]] Hansbo P, Szepessy A. A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations. ''Comput. Methods Appl. Mech. Engrg.'', 84, 175&#8211;192, (1990).
1068
1069
<div id="cite-10"></div>
1070
[[#citeF-10|[10]]]   Cruchaga MA, Oñate E. A generalized streamline finite element approach for the analysis of incompressible flow problems including moving surfaces. ''Comput. Methods in Appl. Mech. Engrg.'', 173, 241&#8211;255, (1999).
1071
1072
<div id="cite-11"></div>
1073
[[#citeF-11|[11]]] Hughes TJR, Franca LP, Hulbert GM. A new finite  element formulation for computational fluid dynamics: VIII. The  Galerkin/least-squares method for advective-diffusive equations. '' Comput. Methods Appl. Mech. Engrg.'', 73, 173&#8211;189, (1989).
1074
1075
<div id="cite-12"></div>
1076
[[#citeF-12|[12]]] Tezduyar TE, Mittal S, Ray SE, Shih R. Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity&#8211;pressure elements. ''Comput. Methods Appl. Mech. Engrg.'', 95, 221&#8211;242,  (1992).
1077
1078
<div id="cite-13"></div>
1079
[[#citeF-13|[13]]]  Douglas J, Russell TF. Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. ''SIAM J. Numer. Anal.'', 19, 871, (1982).
1080
1081
<div id="cite-14"></div>
1082
[[#citeF-14|[14]]]  Pironneau O. On the transport-diffusion algorithm and its applications to the Navier-Stokes equations. ''Numer. Math.'', 38, 309, (1982)..
1083
1084
<div id="cite-15"></div>
1085
[[#citeF-15|[15]]]  Löhner R., Morgan K, Zienkiewicz OC. The solution of non-linear hyperbolic equation systems by the finite element method. ''Int. J. Num. Meth. in Fluids'', 4, 1043, (1984).
1086
1087
<div id="cite-16"></div>
1088
[[#citeF-16|[16]]] Zienkiewicz  OC,  Codina R. A general algorithm for  compressible and incompressible flow. Part I: The split characteristic  based scheme. ''Int. J. Num. Meth. in Fluids'', 20, 869-85, (1995).
1089
1090
<div id="cite-17"></div>
1091
[[#citeF-17|[17]]]  Codina R, Vazquez M, Zienkiewicz OC. A general algorithm for compressible and incompressible flow - Part III. The semi-implicit form. ''Int. J. Num. Meth. in Fluids'', 27, 13&#8211;32, (1998).
1092
1093
<div id="cite-18"></div>
1094
[[#citeF-18|[18]]] Hughes TJR. Multiscale phenomena: Green functions,  subgrid scale models, bubbles and the origins of stabilized methods.  ''Comput. Methods Appl. Mech. Engrg'', 127, 387&#8211;401, (1995).
1095
1096
<div id="cite-19"></div>
1097
[[#citeF-19|[19]]] Brezzi F, Franca LP, Hughes TJR, Russo A. <math>b=\int  g</math>. ''Comput. Methods Appl. Mech. Engrg.'', 145, 329&#8211;339, (1997).
1098
1099
<div id="cite-20"></div>
1100
[[#citeF-20|[20]]]  Codina R. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element method. ''Comput. Methods Appl. Mech. Engrg.'', 190, 1579&#8211;1599, (2000).
1101
1102
<div id="cite-21"></div>
1103
[[#citeF-21|[21]]]  Hauke G. A simple subgrid scale stabilized method for the advection-diffusion-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', 191, 2925&#8211;2948, (2002).
1104
1105
<div id="cite-22"></div>
1106
[[#citeF-22|[22]]] Oñate E. Derivation of stabilized equations for  advective-diffusive transport and fluid flow problems.  ''Comput. Methods Appl. Mech. Engrg.'', 151:1-2, 233&#8211;267, (1998).
1107
1108
<div id="cite-23"></div>
1109
[[#citeF-23|[23]]] Oñate E, Manzan M. Stabilization techniques for finite element analysis of convection diffusion problems. In ''Comput. Anal. of Heat Transfer'',  WIT Press, G. Comini and B. Sunden (Eds.), (2000).
1110
1111
<div id="cite-24"></div>
1112
[[#citeF-24|[24]]] Oñate E, García J, Idelsohn SR. Computation of the stabilization  parameter for the finite element solution of advective-diffusive problems. ''Int. J. Num. Meth. Fluids'', 25, 1385&#8211;1407, (1997).
1113
1114
<div id="cite-25"></div>
1115
[[#citeF-25|[25]]] Oñate E, García J, Idelsohn SR.  An alpha-adaptive  approach for stabilized finite element solution of advective-diffusive  problems with sharp gradients. ''New Advances in Adaptive  Comput. Methods in Mech.'', Elsevier, P. Ladeveze, J.T. Oden (Eds.), (1998).
1116
1117
<div id="cite-26"></div>
1118
[[#citeF-26|[26]]] Oñate E, Manzan M. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. ''Int. J. Num. Meth. Fluids'', 31, 203&#8211;221, (1999).
1119
1120
<div id="cite-27"></div>
1121
[[#citeF-27|[27]]] Oñate E. Possibilities of finite calculus in computational mechanics. ''Int. J. Num. Meth. Engng.'', (2004).
1122
1123
<div id="cite-28"></div>
1124
[[#citeF-28|[28]]]  Oñate E, Taylor RL, Zienkiewicz OC,   Rojek J. A residual correction method based on finite calculus. ''Engineering Computatlions'', 20:5/6, 629&#8211;658, (2003).
1125
1126
<div id="cite-29"></div>
1127
[[#citeF-29|[29]]]  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Comput. Methods Appl. Mech. Engrg.'', 182:1&#8211;2, 355&#8211;370, (2000).
1128
1129
<div id="cite-30"></div>
1130
[[#citeF-30|[30]]]  Oñate E,  García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in ''Comput. Methods Appl. Mech. Engrg.'', 191:6-7, 635-660, (2001).
1131
1132
<div id="cite-31"></div>
1133
[[#citeF-31|[31]]]  García J, Oñate E. An unstructured finite element solver for ship hydrodynamics. ''J. Appl. Mechanics'', 70, 18&#8211;26, (2003).
1134
1135
<div id="cite-32"></div>
1136
[[#citeF-32|[32]]]  Oñate E, García J, Bugeda G, Idelsohn SR. A general stabilized formulation for incompressible fluid flow using finite calculus and the FEM. In ''Towards a New Fluid Dynamics with its Challenges in Aeronautics'', J. Periaux, D. Champion, O. Pironneau and Ph. Thomas (Eds.), CIMNE, Barcelona, (2002).
1137
1138
<div id="cite-33"></div>
1139
[[#citeF-33|[33]]]   Oñate E, García J, Idelsohn SR. Ship Hydrodynamics. ''Encyclopedia of Computational Mechanics'', T. Hughes, R. de Borst and E. Stein (Eds.), J. Wiley, (2004) (to be published).
1140
1141
<div id="cite-34"></div>
1142
[[#citeF-34|[34]]]  Hughes TJR, Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system. ''Comput. Methods Appl. Mech. Engrg.'', 58, 329&#8211;336, (1986b).
1143
1144
<div id="cite-35"></div>
1145
[[#citeF-35|[35]]]  Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation. ''Comput. Methods Appl. Mech. Engrg.'', 110, 325&#8211;342, (1993).
1146
1147
<div id="cite-36"></div>
1148
[[#citeF-36|[36]]]  Tezduyar T. Stabilization parameters and local length scales in SUPG and PSPG formulations. Proceedings of the ''Fith World Congress on Computational Mechanics'', (http://wccm.tuwien.ac.at), Vienna, Austria, July 7&#8211;12 (2002).
1149
1150
<div id="cite-37"></div>
1151
[[#citeF-37|[37]]]  Zienkiewicz OC, Zhu JZ. The Superconvergent patch recovery (SPR) and adaptive finite element refinement. ''Comput. Methods Appl. Mech. Engrg.'',  101, 207&#8211;224, (1992).
1152
1153
<div id="cite-38"></div>
1154
[[#citeF-38|[38]]]  Wiberg NW, Abdulwahab F,  Li XD. (1997). Error estimation  and adaptive procedures based on superconvergent patch recovery. ''Archives Comput. Meth. Engng.'', 4:3, 203&#8211;242, (1997).
1155
1156
<div id="cite-39"></div>
1157
[[#citeF-39|[39]]]  Ilinca F, Hétu JF, Pelletier D. On stabilized finite element formulation for incompressible advective-diffusive transport and fluid flow problems. ''Comp. Methods Appl. Mech. Engrg.'', 188, 235&#8211;257, (2000).
1158
1159
<div id="cite-40"></div>
1160
[[#citeF-40|[40]]]  Chetverushkin BN. Kinetic schemes and their application for simulation of viscous gas dynamics problems. ''European Congres on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2000)'', Barcelona 11&#8211;14 September, E. Oñate ''et al.'' (Eds.), (2000).
1161
1162
<div id="cite-41"></div>
1163
[[#citeF-41|[41]]]  Oñate E, Rojek J, Taylor RL, Zienkiewicz OC. Linear triangles and tetrahedra for incompressible problems using a finite calculus formulation. In ''European Conference on Computational Mechanics (ECCM2001)'', Cracow, Poland, 26&#8211;29 June, (2001).
1164
1165
<div id="cite-42"></div>
1166
[[#citeF-42|[42]]] Oñate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. To appear in ''Int. J. Num. Meth. Engng.'', (2004).
1167
1168
<div id="cite-43"></div>
1169
[[#citeF-43|[43]]]  Oñate E, Zarate F and Idelsohn SR. Stabilized FEM for convection-diffusion problems using a FIC approach with gradient-based stabilization parameters. Submitted to ''Comput. Meth. Appl. Mech. Eng.'', (2004).
1170
1171
<div id="cite-44"></div>
1172
[[#citeF-44|[44]]]  Oñate E and Flores F. Refined finite calculus method for finite element analysis of incompressible flows. Submitted to ''Comput. Meth. Appl. Mech. Eng.'', (2004).
1173
1174
<div id="cite-45"></div>
1175
[[#citeF-45|[45]]] Oñate E and Rojek, J. FIC formulation for finite element analysis of incompressible solids. Submitted to ''Int. J. Num. Meth. Eng.'', (2004).
1176

Return to Onate 2004a.

Back to Top