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
==STABILIZED SOLUTION  OF  THE MULTIDIMENSIONAL ADVECTION-DIFFUSION-ABSORPTION EQUATION USING  LINEAR FINITE ELEMENTS==
2
3
'''Eugenio Oñate, Juan Miquel  and Francisco Zárate'''
4
5
==Abstract==
6
7
A stabilized finite element method (FEM) for the multidimensional steady state advection-diffusion-absorption equation is presented. The stabilized formulation is based on the modified governing differential equations derived via the Finite Calculus (FIC) method. For 1D problems the stabilization terms act as a nonlinear additional diffusion governed by a single stabilization parameter. It is shown that for multidimensional problems an orthotropic stabilizing diffusion must be added along the principal directions of  curvature of the solution. A simple iterative algorithm yielding a stable and accurate solution for all the range of physical parameters and boundary conditions is described. Numerical results for 1D and 2D problems with sharp gradients are presented showing the effectiveness and accuracy of the new stabilized formulation.
8
9
==1 INTRODUCTION==
10
11
Considerable effort has been spent in recent years to derive finite element methods (FEM) <span id='citeF-1'></span>[[#cite-1|1]] for the solution of the advection-diffusion-reaction equation. In this work we will focus on the so called ''exponential regime'' originated by large absorptive (dissipative) reaction terms. Here the solutions are of the form of real exponential functions.  Numerical schemes  find difficulties to approximating the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damköhler numbers. Non physical oscilaltory solution are found with the standard Galerkin FEM unless some stabilization procedure is used.
12
13
Stabilized methods to tackle this problem have been based on  streamline-upwind/Petrov-Galerkin (SUPG)  <span id='citeF-2'></span>[[#cite-2|2]], Galerkin/least-squares <span id='citeF-5'></span>[[#cite-5|5]], Subgrid Scale (SGS) <span id='citeF-5'></span>[[#cite-5|5]] and Residual Free Bubbles <span id='citeF-14'></span>[[#cite-14|14]] finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection-diffusion and the diffusion-reaction equations (Vol. 3 in <span id='citeF-1'></span>[[#cite-1|1]] and <span id='citeF-8'></span>[[#cite-8|8]]), this is not the case for the diffusion-advection-reaction equation. Here, in general,  ''two stabilization parameters'' are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions <span id='citeF-4'></span>[[#cite-4|4]]. As reported in <span id='citeF-12'></span>[[#cite-12|12]] the SUPG, GLS and SGS methods with a single stabilization parameter fail to obtain a stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution.
14
15
Oñate ''et al.'' [18] have recently presented a stabilized FEM for the advection-diffusion absorption equation based on the use of a single stabilization parameter which has the form of a diffusion term. In [18] the formulation is detailed for 1D problems and only a brief introduction to the multidimensional case is given. This paper extends the ideas presented in [18] and provides evidence of the effectiveness and accuracy of the new formulation to deal with multidimensional advection-diffusion-absorption  problems  with sharp gradients.
16
17
The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the ''Finite Calculus'' (FIC) method [19&#8211;20]. The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. Although the FIC&#8211;FEM formulation here presented is general, we will restrict its application in this work to linear finite element approximations only.
18
19
The Galerkin FIC-FEM formulation  described here introduces naturally an additional nonlinear dissipation term into the discretized equations which is governed by a ''single stabilization parameter''. In the absence of the absorption term the formulation simplifies to the standard Petrov-Galerkin approach for the advection-diffusion problem For the diffusion-absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and Oñate using a variational FIC approach  <span id='citeF-15'></span>[[#cite-15|15]]. The general nonlinear form of the stabilization parameter  is a function of the signs of the solution and   its first and second derivatives. This introduces a non-linearity in the solution scheme and a simple iterative algorithm  is described. A simpler constant expression of the stabilization parameter is also presented.
20
21
Details of the 1D formulation and its extension to deal with multidimensional problems are given. For the multidimensional case Oñate ''et al.'' <span id='citeF-27'></span>[[#cite-27|27]] have recently shown that a general form of the stabilization parameters can be found by writting the FIC equations along the principal curvature directions of the solution. The resulting FIC-FEM formulation is equivalent in this case (for linear elements) to adding a stabilizing diffusion matrix to the standard infinitessimal equation. The stabilizing diffusion matrix depends on the signs of the solution and its derivatives and on the velocities along the principal curvature directions of the solution. This introduces a nonlinearity in the solution process. We present a simple iterative scheme based in assuming that  the main principal curvature direction at each  point is coincident with the gradient vector direction. In the last part of the paper we present a collection of 1D and 2D examples showing the effectiveness and accuracy of the new FIC-FEM formulation for different values of the advective and absorptive terms.
22
23
==2 FIC FORMULATION OF THE 1D STATIONARY ADVECTION-DIFFUSION-ABSORPTION EQUATION==
24
25
The governing equation for the 1D stationary advection-diffusion-absorption problem can be written in the FIC formulation as
26
27
{| class="formulaSCP" style="width: 100%; text-align: left;" 
28
|-
29
| 
30
{| style="text-align: left; margin:auto;" 
31
|-
32
| style="text-align: center;" | <math>r - \underline{{h\over 2} {dr\over dx}}{h\over 2} {dr\over dx}=0\quad \hbox{in } x \in (0,L) </math>
33
|}
34
| style="width: 5px;text-align: right;" | (1)
35
|}
36
37
{| class="formulaSCP" style="width: 100%; text-align: left;" 
38
|-
39
| 
40
{| style="text-align: left; margin:auto;" 
41
|-
42
| style="text-align: center;" | <math>-u\phi + k {d\phi \over dx} +q^p - \underline{{h\over 2} r}{h\over 2} r=0\quad \hbox{on }\Gamma _q  </math>
43
|}
44
| style="width: 5px;text-align: right;" | (2)
45
|}
46
47
{| class="formulaSCP" style="width: 100%; text-align: left;" 
48
|-
49
| 
50
{| style="text-align: left; margin:auto;" 
51
|-
52
| style="text-align: center;" | <math>\phi -\phi ^p =0 \quad \hbox{on }\Gamma _\phi  </math>
53
|}
54
| style="width: 5px;text-align: right;" | (3)
55
|}
56
57
where
58
59
{| class="formulaSCP" style="width: 100%; text-align: left;" 
60
|-
61
| 
62
{| style="text-align: left; margin:auto;" 
63
|-
64
| style="text-align: center;" | <math>r =:-u {d\phi \over dx} +{d\over dx} \left(k{d\phi \over dx}\right)- s\phi + Q  </math>
65
|}
66
| style="width: 5px;text-align: right;" | (4)
67
|}
68
69
In above equations <math display="inline">\phi </math> is the state variable, <math display="inline">x \in [0,L]</math> is the problem domain, <math display="inline">L</math> is the domain length, <math display="inline">u</math> is the velocity field, <math display="inline">k\ge 0</math> is the diffusion, <math display="inline">s\ge 0</math> is the absorption, dissipation or destruction source parameter, <math display="inline">Q</math> is a constant source term, <math display="inline">q^p</math> and <math display="inline">\phi ^p</math> are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries <math display="inline">\Gamma _q</math> and <math display="inline">\Gamma _\phi </math>, respectively and <math display="inline">h</math> is a ''characteristic length'' which plays the role of a stabilization parameter. In the 1D problem <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math> consist of four combinations at the end points of the problem domain.
70
71
Eqs.(1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitessimal theory [19,20]. The underlined stabilizing terms in Eqs.(1) and (2) emanate from these higher order expansions. Note that as the characteristic length parameter <math display="inline">h</math> tends to zero the FIC differential equations gradually recover the standard infinitessimal form.
72
73
Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [19&#8211;30,37].
74
75
==3 FINITE ELEMENT FORMULATION==
76
77
We will construct a standard finite element discretization <math display="inline">\left\{l^e\right\}</math> of the 1D analysis domain length <math display="inline">L</math> with index <math display="inline">e</math> ranging from 1 to the number of elements <math display="inline">N</math> <span id='citeF-1'></span>[[#cite-1|1]]. The state variable <math display="inline">\phi </math> is approximated by <math display="inline">\bar \phi </math> over the analysis domain. The approximated variable <math display="inline">\bar \phi </math> is interpolated within each element with <math display="inline">n</math> nodes in the standard manner, i.e.
78
79
{| class="formulaSCP" style="width: 100%; text-align: left;" 
80
|-
81
| 
82
{| style="text-align: left; margin:auto;" 
83
|-
84
| style="text-align: center;" | <math>\phi \simeq \bar \phi \quad \hbox{for}\quad x \in [0,L]</math>
85
|}
86
| style="width: 5px;text-align: right;" |  (5a)
87
|}
88
89
with
90
91
{| class="formulaSCP" style="width: 100%; text-align: left;" 
92
|-
93
| 
94
{| style="text-align: left; margin:auto;" 
95
|-
96
| style="text-align: center;" | <math>\bar \phi =\sum \limits _{i=1}^n N_i \phi _i</math>
97
|}
98
| style="width: 5px;text-align: right;" |  (5b)
99
|}
100
101
where <math display="inline">N_i</math> are the element shape functions and <math display="inline">\phi _i</math> are nodal values of the approximate function <math display="inline">\bar \phi </math>. Substituting Eq.(5a) into Eqs.(1) and (2) gives
102
103
{| class="formulaSCP" style="width: 100%; text-align: left;" 
104
|-
105
| 
106
{| style="text-align: left; margin:auto;" 
107
|-
108
| style="text-align: center;" | <math>\bar r - {h\over 2} {d\bar r\over dx} =r_\Omega \quad \hbox{in } x\in (0,L) </math>
109
|}
110
|}
111
112
{| class="formulaSCP" style="width: 100%; text-align: left;" 
113
|-
114
| 
115
{| style="text-align: left; margin:auto;" 
116
|-
117
| style="text-align: center;" | <math>-u\bar \phi + k {d\bar \phi \over dx} +q^p - {h\over 2} \bar r=r_q\quad \hbox{on }\Gamma _q  </math>
118
|}
119
| style="width: 5px;text-align: right;" | (6)
120
|}
121
122
{| class="formulaSCP" style="width: 100%; text-align: left;" 
123
|-
124
| 
125
{| style="text-align: left; margin:auto;" 
126
|-
127
| style="text-align: center;" | <math>\bar \phi -\phi ^p =r_\phi \quad \hbox{on }\Gamma _\phi  </math>
128
|}
129
| style="width: 5px;text-align: right;" | (7)
130
|}
131
132
where <math display="inline">\bar r =r(\bar \phi )</math> and <math display="inline">r_\Omega , r_q</math> and <math display="inline">r_\phi </math> are the residuals of the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries <math display="inline">\Gamma _q</math> and <math display="inline">\Gamma _\phi </math>, respectively.
133
134
The weighted residual form of Eqs.(6)&#8211;(8) is written as
135
136
{| class="formulaSCP" style="width: 100%; text-align: left;" 
137
|-
138
| 
139
{| style="text-align: left; margin:auto;" 
140
|-
141
| style="text-align: center;" | <math>\int _L W_i \left(\bar r - {h\over 2} {d\bar r\over dx}\right)dx + \left[\hat W_i \left(-u \bar \phi + k {d\bar \phi \over dx} +q_p - {h\over 2} \bar r \right)\right]_{\Gamma _q} =0 </math>
142
|}
143
| style="width: 5px;text-align: right;" | (8)
144
|}
145
146
where <math display="inline">W_i(x)</math> and <math display="inline">\hat W_i</math> are test functions satisfying <math display="inline">W_i =\hat W_i =0</math> on <math display="inline">\Gamma _\phi </math>.
147
148
Assuming smooth enough solutions and integrating by parts the term involving <math display="inline">h</math> in the first integral gives for <math display="inline">\hat W_i=-W_i</math>
149
150
{| class="formulaSCP" style="width: 100%; text-align: left;" 
151
|-
152
| 
153
{| style="text-align: left; margin:auto;" 
154
|-
155
| style="text-align: center;" | <math>\int _L W_i \bar r d x - \left[W_i \left(-u \bar \phi + k {d\bar \phi \over dx} +q^p \right)\right]_{\Gamma _q} +\sum \limits _e \int _{l^e} {h\over 2} {dW_i\over dx} \bar r dx =0</math>
156
|}
157
| style="width: 5px;text-align: right;" | (9)
158
|}
159
160
The third term in Eq.(10) is computed as the sum of the integrals over the element interiors, therefore allowing for the space derivatives of <math display="inline">\bar r</math> to be discontinuous. Also in  Eq.(10) <math display="inline">h</math> has been assumed to be constant within each element, (i.e. <math display="inline">\displaystyle{dh\over dx}=0</math> within <math display="inline">l^e</math>).
161
162
The  weak form is obtained by integrating by parts the advective and diffusive terms within <math display="inline">\bar r</math> in the first integral of Eq.(10). This gives
163
164
{| class="formulaSCP" style="width: 100%; text-align: left;" 
165
|-
166
| 
167
{| style="text-align: left; margin:auto;" 
168
|-
169
| style="text-align: center;" | <math>\int _L \left[u {dW_i\over dx} \bar \phi - {dW_i\over dx}k {d\bar \phi \over dx}- W_i s \bar \phi + W_i Q\right]dx -  [W_i q^p]_{\Gamma _q} -  \sum \limits _e \int _{l^e} \left(\beta k {dW_i\over dx} {d\bar \phi \over dx} -{h\over 2} {dW_i\over dx}Q\right)dx =0 </math>
170
|}
171
| style="width: 5px;text-align: right;" | (10)
172
|}
173
174
with
175
176
{| class="formulaSCP" style="width: 100%; text-align: left;" 
177
|-
178
| 
179
{| style="text-align: left; margin:auto;" 
180
|-
181
| style="text-align: center;" | <math>\beta = \left[{s\bar \phi \over 2k\bar \phi '}+{u\over 2k}-{\bar \phi ''\over 2\bar \phi '} \right]h </math>
182
|}
183
| style="width: 5px;text-align: right;" | (11)
184
|}
185
186
where a prime denotes the derivative with respect to the space coordinate.
187
188
Wee see clearly that the last term of Eq.(11) introduces within each element an additional diffusion of value <math display="inline">\beta k</math>.
189
190
Substituting expression (5b) into (11) and choosing a Galerkin method with <math display="inline">W_i =N_i</math> within each element gives the discrete system of FE equations written in the standard matrix form as
191
192
{| class="formulaSCP" style="width: 100%; text-align: left;" 
193
|-
194
| 
195
{| style="text-align: left; margin:auto;" 
196
|-
197
| style="text-align: center;" | <math>{K}\bar{\boldsymbol \phi } ={f} </math>
198
|}
199
| style="width: 5px;text-align: right;" | (12)
200
|}
201
202
where <math display="inline">\bar{\boldsymbol \phi }</math> is the vector of nodal unknowns and the element contributions to matrix '''K''' and vector <math display="inline">f</math> are
203
204
{| class="formulaSCP" style="width: 100%; text-align: left;" 
205
|-
206
| 
207
{| style="text-align: left; margin:auto;" 
208
|-
209
| style="text-align: right;" | <math>K_{ij}^e\!\! </math>
210
| style="text-align: center;" | <math>=</math>
211
| <math>\!\! \int _{l^e} \left(-u {dN_i\over dx} N_j + k(1+\beta ) {dN_i\over dx}{dN_j\over dx}+ sN_i N_j\right)dx</math>
212
| style="width: 5px;text-align: right;" | (13)
213
|-
214
| style="text-align: right;" | <math> f_i^e\!\! </math>
215
| style="text-align: center;" | <math>=</math>
216
| <math> \!\! \int _{l^e} \left(N_i + {h\over 2} {dN_i\over dx} \right)Qdx - (N_i q^p)_{\Gamma _q} </math>
217
| style="width: 5px;text-align: right;" | (14)
218
|}
219
|}
220
221
The amount of balancing diffusion in Eq.(14) clearly depends on the (nonlinear) stabilization parameter <math display="inline">\beta </math>. The element and critical  values of <math display="inline">\beta </math> are deduced in the next section for linear two node elements.
222
223
We note that the integral of the term <math display="inline">\displaystyle{h\over 2} \displaystyle{dN_i\over dx}Q</math> in Eq.(15) vanishes after asssembly when <math display="inline">h </math> and <math display="inline">Q</math> are uniform over a patch of linear elements.
224
225
==4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS==
226
227
The matrix <math display="inline">{K}^e</math> and the vector <math display="inline">{f}^e</math> for two node linear elements are (for constant values of <math display="inline">u,k</math>, <math display="inline">s</math> and <math display="inline">Q</math>)
228
229
{| class="formulaSCP" style="width: 100%; text-align: left;" 
230
|-
231
| 
232
{| style="text-align: left; margin:auto;" 
233
|-
234
| style="text-align: center;" | <math>{K}^e = - {u\over 2} \left[\begin{matrix}-1 & -1\\ 1 & 1\\\end{matrix}\right]+ {k\over l^e} (1+\beta ^e)  \left[\begin{matrix}1 & -1\\ - 1 & 1\\\end{matrix}\right]+ {sl^e \over 6} \left[\begin{matrix}2 & 1\\ 1 & 2\\\end{matrix}\right]</math>
235
|}
236
| style="width: 5px;text-align: right;" |  (16a)
237
|}
238
239
{| class="formulaSCP" style="width: 100%; text-align: left;" 
240
|-
241
| 
242
{| style="text-align: left; margin:auto;" 
243
|-
244
| style="text-align: center;" | <math>{f}^e = {Ql^e\over 2}\left\{\begin{matrix}1 - \displaystyle{h^e\over 2}\\ 1+ \displaystyle{h^e\over 2}\\\end{matrix}\right\}\qquad + \quad \hbox{boundary term}</math>
245
|}
246
| style="width: 5px;text-align: right;" |  (16a)
247
|}
248
249
In Eqs.(16) index <math display="inline">e</math>denotes element values.
250
251
Assuming <math display="inline">Q=0</math>, a typical stencil for elements of equal size <math display="inline">l</math> can be written as
252
253
{| class="formulaSCP" style="width: 100%; text-align: left;" 
254
|-
255
| 
256
{| style="text-align: left; margin:auto;" 
257
|-
258
| style="text-align: center;" | <math>\begin{array}{r}-\gamma (\bar \phi _{i+1} -\bar \phi _{i-1})-(1+\beta ) \bar \phi _{i-1}+2(1+\beta ) \bar \phi _i - (1+\beta ) \bar \phi _{i+1}+\\ + \displaystyle{w\over 6} (\bar \phi _{i-1}+ 4 \bar \phi _i + \bar \phi _{i+1})=0 \end{array} </math>
259
|}
260
| style="width: 5px;text-align: right;" | (17)
261
|}
262
263
where for simplicity a constant value of <math display="inline">\beta </math> across the mesh has been assumed. In Eq.(17) <math display="inline">\gamma ={ul\over 2k}</math> and <math display="inline">w={sl^2\over 4}</math> are the Peclet number and a velocity independent dimensionless number, respectively.
264
265
From Eq.(17) we deduce
266
267
{| class="formulaSCP" style="width: 100%; text-align: left;" 
268
|-
269
| 
270
{| style="text-align: left; margin:auto;" 
271
|-
272
| style="text-align: center;" | <math>\beta = \gamma \left({\bar \phi _{i+1} -\bar \phi _{i-1}\over \bar \phi _{i+1}- 2\bar \phi _i+  \bar \phi _{i-1} }\right)+ {w\over 6} \left({\bar \phi _{i-1} + 4 \bar \phi _i+\bar \phi _{i+1}\over \bar \phi _{i+1}- 2\bar \phi _i+ \bar \phi _{i+1} }\right)-1 </math>
273
|}
274
| style="width: 5px;text-align: right;" | (18)
275
|}
276
277
In the vecinity of a sharp gradient zone we can take
278
279
{| class="formulaSCP" style="width: 100%; text-align: left;" 
280
|-
281
| 
282
{| style="text-align: left; margin:auto;" 
283
|-
284
| style="text-align: center;" | <math>\begin{array}{l}\bar \phi _{i+1} -\bar \phi _{i-1}\simeq \bar \phi _{max} S_1\\ \bar \phi _{i+1} -2\bar \phi _i+  \bar \phi _{i-1}=\bar \phi _{max} S_2\\ \bar \phi _i+4 \bar \phi _i + \bar \phi _{i+1}= \bar \phi _{i+1}S_0 \end{array} </math>
285
|}
286
| style="width: 5px;text-align: right;" | (19)
287
|}
288
289
where <math display="inline">\bar \phi _{max}</math> is the maximum value of the approximate function <math display="inline">\bar \phi </math> in the patch of elements adjacent to the sharp gradient zone and
290
291
{| class="formulaSCP" style="width: 100%; text-align: left;" 
292
|-
293
| 
294
{| style="text-align: left; margin:auto;" 
295
|-
296
| style="text-align: center;" | <math>S_0= \hbox{sign } (\bar \phi ) ,\quad S_1 = \hbox{sign } \left({d\bar \phi \over dx}\right),\quad S_2 = \hbox{sign } \left({d^2 \bar \phi \over dx^2}\right) </math>
297
|}
298
| style="width: 5px;text-align: right;" | (20)
299
|}
300
301
where sign <math display="inline">\bar{(\cdot )}</math> denotes the sign of the magnitude within the brackets computed at the  patch mid point.
302
303
Substituting Eq.(19) into (18) leads to the following expression of the stabilization parameter
304
305
{| class="formulaSCP" style="width: 100%; text-align: left;" 
306
|-
307
| 
308
{| style="text-align: left; margin:auto;" 
309
|-
310
| style="text-align: center;" | <math>\beta = \left[\left({S_0\over S_2}\right){w\over 6} +\left({S_1\over S_2}\right)  \gamma -1\right] </math>
311
|}
312
| style="width: 5px;text-align: right;" | (21)
313
|}
314
315
The ''element stabilization parameter'' <math display="inline">\beta ^e</math> is now defined as
316
317
{| class="formulaSCP" style="width: 100%; text-align: left;" 
318
|-
319
| 
320
{| style="text-align: left; margin:auto;" 
321
|-
322
| style="text-align: center;" | <math>\begin{array}{l}\beta ^e =\beta \quad \hbox{for } \beta >0\\ \beta ^e =0 \quad \hbox{for } \beta \le 0 \end{array} </math>
323
|}
324
| style="width: 5px;text-align: right;" | (22)
325
|}
326
327
where <math display="inline">\beta </math> is given by Eq.(21) and the signs <math display="inline">S_0</math>, <math display="inline">S_1</math> and <math display="inline">S_2</math> are computed now at the element mid-point.
328
329
It is clear from above that the computation of the stabilization parameter <math display="inline">\beta ^e</math> requires the knowledge of the sign of the numerical solution <math display="inline">\bar \phi </math> and that of the  first and second derivatives of <math display="inline">\bar \phi </math> within each element. This necessarily leads to an iterative scheme. A simple algorithm which provides a stabilized and accurate solution in just two steps is presented below.
330
331
===4.1 Critical stabilization parameter and unstability conditions===
332
333
The following constant value of <math display="inline">\beta </math> over the mesh ensures a stabilized solution for all ranges of <math display="inline">\gamma </math> and <math display="inline">w</math>
334
335
{| class="formulaSCP" style="width: 100%; text-align: left;" 
336
|-
337
| 
338
{| style="text-align: left; margin:auto;" 
339
|-
340
| style="text-align: center;" | <math>\displaystyle{\beta \le \beta _c = {w\over 6} +\vert \gamma \vert -1} </math>
341
|}
342
| style="width: 5px;text-align: right;" | (23)
343
|}
344
345
where <math display="inline">\beta _c </math> is the ''critical stabilization parameter''. Note that <math display="inline">\beta _c</math> corresponds to the maximum value of <math display="inline">\beta </math> in Eq.(21) for <math display="inline">{S_0\over S_2}={S_1\over S_2}=1</math>. A mathematical proof of Eq.(23) is given in [18].
346
347
Clearly the value of <math display="inline">\beta _c</math> of Eq.(23) is meaningful only if <math display="inline">\beta _c >0</math> and this can be taken as an indicator of an unstable solution. Conversely, a value of <math display="inline">\beta _c \le 0</math> indicates that no stabilization is needed.
348
349
===4.2 Iterative solution scheme===
350
351
The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution.
352
353
step step
354
355
Compute a first stabilized solution <math display="inline">\bar{\boldsymbol \phi }^1</math> using the critical value <math display="inline">\beta ^e = \beta _c</math> given by Eq.(23). This ensures a stabilized, although sometimes slightly overdiffusive, solution.
356
357
==~&nbsp;==
358
359
step
360
361
Compute the signs of the first and second derivatives of <math display="inline">\bar{\boldsymbol \phi }^1</math> within each element. The second derivative field is obtained as follows
362
363
{| class="formulaSCP" style="width: 100%; text-align: left;" 
364
|-
365
| 
366
{| style="text-align: left; margin:auto;" 
367
|-
368
| style="text-align: center;" | <math>\left({d^2 \bar \phi ^1 \over dx^2}\right)^e= {1\over l^e} \left[\left({d\hat \phi ^1 \over dx}\right)^e_2 - \left({d\hat \phi ^1 \over dx}\right)^e_1\right] </math>
369
|}
370
| style="width: 5px;text-align: right;" | (24)
371
|}
372
373
where <math display="inline">({\hat \cdot })_i^e</math> denotes averaged values of the first derivative at node <math display="inline">i</math> of element <math display="inline">e</math>. At the boundary nodes the constant value of the derivative of <math display="inline">\bar \phi </math> within the element is taken in Eq.(24); i.e. <math display="inline">(\hat{\cdot })_i^e = \left({d\bar \phi \over dx}\right)^{(e)} = {\bar \phi _2 - \bar \phi _1 \over l^e}</math>.
374
375
Compute the enhanced stabilized solution <math display="inline">{\boldsymbol \phi }^2</math> using the element value of <math display="inline">\beta ^e</math> given by Eq.(22).
376
377
In all the 1D examples solved the above two steps have sufficed to obtain a converged stabilized and accurate solution [18]. The reason of this is that  the signs of the first and second derivative fields typically do not change any further after the second step over the elements adjacent to high gradient zones.
378
379
==5 EXTENSION TO MULTI-DIMENSIONAL PROBLEMS==
380
381
Consider the general steady-state advection-diffusion-reaction equation written for the zero constant source case (<math display="inline">Q=0</math>) as
382
383
<span id="eq-25"></span>
384
{| class="formulaSCP" style="width: 100%; text-align: left;" 
385
|-
386
| 
387
{| style="text-align: left; margin:auto;" 
388
|-
389
| style="text-align: center;" | <math>r(\phi ): =- {u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T {D}{\boldsymbol \nabla }\phi - s\phi =0 </math>
390
|}
391
| style="width: 5px;text-align: right;" | (25)
392
|}
393
394
For 2D problems
395
396
<span id="eq-26"></span>
397
{| class="formulaSCP" style="width: 100%; text-align: left;" 
398
|-
399
| 
400
{| style="text-align: left; margin:auto;" 
401
|-
402
| style="text-align: center;" | <math>{u}=[u,v]^T\quad ,\quad {\boldsymbol \nabla }=\left[{\partial  \over \partial x},{\partial  \over \partial y}\right]^T\quad ,\quad {D} =k \left[\begin{matrix}1 &0\\ 0&1\\\end{matrix}\right] </math>
403
|}
404
| style="width: 5px;text-align: right;" | (26)
405
|}
406
407
are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix.
408
409
The FIC form of Eq.([[#eq-25|25]]) is written as
410
411
<span id="eq-27"></span>
412
{| class="formulaSCP" style="width: 100%; text-align: left;" 
413
|-
414
| 
415
{| style="text-align: left; margin:auto;" 
416
|-
417
| style="text-align: center;" | <math>r - \underline{{1\over 2} {h}^T {\boldsymbol \nabla }r}{1\over 2} {h}^T {\boldsymbol \nabla }r=0 </math>
418
|}
419
| style="width: 5px;text-align: right;" | (27)
420
|}
421
422
where <math display="inline">r</math> is the original infinitessimal differential equation as defined in Eq.([[#eq-25|25]]).
423
424
The Dirichlet and boundary conditions of the FIC formulation are
425
426
<span id="eq-28"></span>
427
{| class="formulaSCP" style="width: 100%; text-align: left;" 
428
|-
429
| 
430
{| style="text-align: left; margin:auto;" 
431
|-
432
| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on}\quad \Gamma _\phi   </math>
433
|}
434
| style="width: 5px;text-align: right;" | (28)
435
|}
436
437
<span id="eq-29"></span>
438
{| class="formulaSCP" style="width: 100%; text-align: left;" 
439
|-
440
| 
441
{| style="text-align: left; margin:auto;" 
442
|-
443
| style="text-align: center;" | <math>- {u}^T {n}\phi + {n}^T {D} {\boldsymbol \nabla } \phi + q^p - \underline{{1\over 2} {h}^T {n}r}{1\over 2} {h}^T {n}r=0 \quad \hbox{on}\quad \Gamma _q </math>
444
|}
445
| style="width: 5px;text-align: right;" | (29)
446
|}
447
448
where <math display="inline">n</math> is the normal vector to the boundary where the normal flux is prescribed. As usual index <math display="inline">p</math> denotes the prescribed values.
449
450
In Eqs.([[#eq-27|27]]) and ([[#eq-29|29]]) <math display="inline">{h}=[h_x,h_y]^T</math> is the characteristic vector of the 2D FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs.([[#eq-27|27]]) and ([[#eq-29|29]]) introduce the necessary stability in the numerical solution [19,20,21].
451
452
If vector '''h'''  is taken to be parallel to the velocity '''u''' the standard SUPG method is recovered [18&#8211;23]. The more general form of '''h''' allows to obtain stabilized finite element solutions in the presence of strong gradients of <math display="inline">\phi </math> near the boundaries (boundary layers) and within the analysis domain (internal layers) <span id='citeF-27'></span>[[#cite-27|27]].  The FIC formulation therefore reproduces the best features of the so called shock-capturing or transverse-dissipation schemes <span id='citeF-2'></span>[[#cite-2|2]].
453
454
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
455
|-
456
|[[Image:draft_Samper_447243531-Figure1.png|600px|Global axes (x,y) and principal curvature axes (ξ,η)]]
457
|- style="text-align: center; font-size: 75%;"
458
| colspan="1" | Global axes (<math>x,y</math>) and principal curvature axes (<math>\xi ,\eta </math>)
459
|}
460
461
Let us write down the FIC balance equation  in the principal curvature axes of the solution <math display="inline">\xi ,\eta </math> (Figure 1). The FIC balance equation is
462
463
{| class="formulaSCP" style="width: 100%; text-align: left;" 
464
|-
465
| 
466
{| style="text-align: left; margin:auto;" 
467
|-
468
| style="text-align: right;" | 
469
| style="text-align: center;" | 
470
| <math>-u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+k \left({\partial ^2\phi \over \partial \xi ^2}+ {\partial ^2\phi \over \partial \eta ^2}\right)-s\phi - {h_\xi \over 2} {\partial  \over \partial \xi }\left[-u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+ k \left({\partial ^2\phi \over \partial \xi ^2}+ {\partial ^2\phi \over \partial \eta ^2}\right)-s\phi \right]</math>
471
|-
472
| style="text-align: right;" | 
473
| style="text-align: center;" | 
474
| <math> - {h_\eta \over 2} {\partial  \over \partial \eta } \left[-u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+ k \left({\partial ^2\phi \over \partial \xi ^2}+ {\partial ^2\phi \over \partial \eta ^2}\right)-s\phi \right]=0 </math>
475
|}
476
| style="width: 5px;text-align: right;" | (30)
477
|}
478
479
where <math display="inline">u_\xi , u_\eta </math> are the velocities along the principal axes of curvature <math display="inline">\xi </math> and <math display="inline">\eta </math>, respectively.
480
481
As <math display="inline">\xi </math> and <math display="inline">\eta </math> are the principal curvature axes of the solution then
482
483
{| class="formulaSCP" style="width: 100%; text-align: left;" 
484
|-
485
| 
486
{| style="text-align: left; margin:auto;" 
487
|-
488
| style="text-align: center;" | <math>{\partial ^2\phi \over \partial \xi \partial \eta }= {\partial ^2\phi \over \partial \eta \partial \xi }=0 </math>
489
|}
490
| style="width: 5px;text-align: right;" | (31)
491
|}
492
493
Introducing this simplification into Eq.(30) we can rewrite this equation as
494
495
{| class="formulaSCP" style="width: 100%; text-align: left;" 
496
|-
497
| 
498
{| style="text-align: left; margin:auto;" 
499
|-
500
| style="text-align: center;" | <math>\begin{array}{l} -u_\xi \displaystyle {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+\left(k + {u_\xi h_\xi \over 2}+ {sh_\xi \over 2} {\partial \phi \over \partial \xi } \left({\partial ^2\phi \over \partial \xi ^2} \right)^{-1}  \right){\partial ^2\phi \over \partial \xi ^2} +\\ \displaystyle + \left(k + {u_\eta h_\eta \over 2}+ {sh_\eta \over 2} {\partial \phi \over \partial \eta } \left({\partial ^2\phi \over \partial \eta ^2} \right)^{-1}\right){\partial ^2\phi \over \partial \eta ^2} - s\phi - k \left({h_\xi \over 2}{\partial ^3 \phi \over \partial \xi ^3}+  {h_\eta \over 2} {\partial ^3 \phi \over \partial \eta ^3}\right)=0 \end{array}</math>
501
|}
502
| style="width: 5px;text-align: right;" |  (32a)
503
|}
504
505
or
506
507
{| class="formulaSCP" style="width: 100%; text-align: left;" 
508
|-
509
| 
510
{| style="text-align: left; margin:auto;" 
511
|-
512
| style="text-align: center;" | <math> -u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+ k(1+\beta _\xi ) {\partial ^2\phi \over \partial \xi ^2} + k (1+\beta _\eta ) {\partial ^2\phi \over \partial \eta ^2} - s\phi -k \left({h_\xi \over 2}{\partial ^3 \phi \over \partial \xi ^3}+  {h_\eta \over 2} {\partial ^3 \phi \over \partial \eta ^3}\right)=0</math>
513
|}
514
| style="width: 5px;text-align: right;" |  (32b)
515
|}
516
517
We can see clearly from Eq.(33) that the FIC governing equations introduce orthotropic diffusion parameters  of values <math display="inline">{\beta _\xi k }</math> and <math display="inline"> { \beta _\eta k}</math> along the <math display="inline">\xi </math> and <math display="inline">\eta </math> axes, respectively with
518
519
{| class="formulaSCP" style="width: 100%; text-align: left;" 
520
|-
521
| 
522
{| style="text-align: left; margin:auto;" 
523
|-
524
| style="text-align: center;" | <math>\beta _\xi = {u_\xi h_\xi \over 2k} + {sh_\xi \over 2k} {\partial \phi  \over \partial \xi } \left({\partial ^2\phi \over \partial \xi ^2} \right)^{-1},\quad \beta _\eta ={u_\xi h_\xi \over 2k} + {sh_\eta \over 2k} {\partial \phi  \over \partial \eta } \left({\partial ^2\phi \over \partial \eta ^2} \right)^{-1} </math>
525
|}
526
| style="width: 5px;text-align: right;" |  (33)
527
|}
528
529
Also note that the last term of Eq.(32b) will vanish after discretization for linear elements.
530
531
Eq.(32b) can be rewritten in matrix form (neglecting the last term) as
532
533
{| class="formulaSCP" style="width: 100%; text-align: left;" 
534
|-
535
| 
536
{| style="text-align: left; margin:auto;" 
537
|-
538
| style="text-align: center;" | <math>- {u}^{\prime T} {\boldsymbol \nabla }^\prime \phi + {\boldsymbol \nabla }^{\prime T} ({D}+\bar {D}^\prime ) {\boldsymbol \nabla }^{\prime }\phi - s\phi =0 </math>
539
|}
540
|}
541
542
where <math display="inline">{u}^\prime =[u_\xi , u_\eta ]^T</math>, <math display="inline">{\boldsymbol \nabla }^\prime = \left[{\partial \over \partial \xi }, {\partial \over \partial \eta }\right]^T</math>, <math display="inline">D</math> is the “physical” isotropic diffusion matrix and <math display="inline">\bar {D}'</math> is the balancing diffusion matrix in the local axes <math display="inline">\xi </math> and <math display="inline">\eta </math>. The form of this matrix is
543
544
{| class="formulaSCP" style="width: 100%; text-align: left;" 
545
|-
546
| 
547
{| style="text-align: left; margin:auto;" 
548
|-
549
| style="text-align: center;" | <math>\bar {D}^\prime = \left[\begin{array}{cc}\beta _\xi k& 0\\ 0 &\beta _\eta k \end{array}\right] </math>
550
|}
551
| style="width: 5px;text-align: right;" | (34)
552
|}
553
554
The velocities  along the principal curvature axes <math display="inline">u_\xi </math> and <math display="inline">u_\eta </math> can be obtained by projecting the cartesian velocities into the principal curvature axes <math display="inline">\xi </math> and <math display="inline">\eta </math> as
555
556
{| class="formulaSCP" style="width: 100%; text-align: left;" 
557
|-
558
| 
559
{| style="text-align: left; margin:auto;" 
560
|-
561
| style="text-align: center;" | <math>{u}'=\left\{\begin{array}{c}u_\xi \\ u_\eta \end{array}\right\}= {T} {u}  \quad \hbox{with}\quad {T}= \left[\begin{array}{cc}c_\alpha & s_\alpha \\ - s_\alpha & c_\alpha \end{array}\right]\quad ,\quad {u}=\left\{\begin{array}{c}u\\v\end{array}\right\} </math>
562
|}
563
| style="width: 5px;text-align: right;" | (35)
564
|}
565
566
where <math display="inline">c_\alpha =\cos \alpha </math>, <math display="inline">s_\alpha =\sin \alpha </math> and <math display="inline">\alpha </math> is the angle which the <math display="inline">\xi </math> axis forms with the <math display="inline">x</math> axis (Figure 1). Note that as the solution is continuous the principal curvature directions <math display="inline">\xi </math> and <math display="inline">\eta </math> are orthogonal.
567
568
The values of <math display="inline">\beta _\xi </math> and <math display="inline">\beta _\eta </math> are  computed by considering the solution of two uncoupled 1D problems along the <math display="inline">\xi </math> and <math display="inline">\eta </math> directions. This gives from Eqs.(21) and (22)
569
570
{| class="formulaSCP" style="width: 100%; text-align: left;" 
571
|-
572
| 
573
{| style="text-align: left; margin:auto;" 
574
|-
575
| style="text-align: center;" | <math>\beta _\xi = \left[\left({S_0\over S_{\xi _2}}\right){w_\xi \over 6} +  \left({S_{\xi _1}\over S_{\xi _2}}\right)\gamma _\xi -1\right]\quad ,\quad \gamma _\xi = {u_\xi l_\xi \over 2k} \quad ,\quad  w_\xi = {sl^2_\xi \over k}</math>
576
|}
577
| style="width: 5px;text-align: right;" |  (37a)
578
|}
579
580
{| class="formulaSCP" style="width: 100%; text-align: left;" 
581
|-
582
| 
583
{| style="text-align: left; margin:auto;" 
584
|-
585
| style="text-align: center;" | <math>\beta _\eta = \left[\left({S_0\over S_{\eta _2}}\right){w_\eta \over 6} +  \left({S_{\eta _1}\over S_{\eta _2}}\right)\gamma _\eta -1\right]\quad ,\quad \gamma _\eta = {u_\eta l_\eta \over 2k} \quad ,\quad  w_\eta = {sl^2_\eta \over k} </math>
586
|}
587
| style="width: 5px;text-align: right;" |  (37b)
588
|}
589
590
where
591
592
{| class="formulaSCP" style="width: 100%; text-align: left;" 
593
|-
594
| 
595
{| style="text-align: left; margin:auto;" 
596
|-
597
| style="text-align: center;" | <math>\begin{array}{l}S_0 =\hbox{sign }(\bar \phi ) \quad ,\quad S_{\xi _1}= \hbox{sign } \left(\displaystyle{\partial \bar \phi \over \partial \xi }\right)\quad ,\quad  S_{\xi _2}= \hbox{sign } \left(\displaystyle{\partial ^2 \bar \phi \over \partial \xi ^2}\right)\\ S_{\eta _1}= \hbox{sign } \left(\displaystyle{\partial \bar \phi \over \partial \eta }\right)\quad ,\quad  S_{\eta _2}= \hbox{sign } \left(\displaystyle{\partial ^2 \bar \phi \over \partial \eta ^2}\right) \end{array} </math>
598
|}
599
| style="width: 5px;text-align: right;" | (38)
600
|}
601
602
and <math display="inline">\bar \phi </math> is as usual the approximate solution.
603
604
The lengths <math display="inline">l_\xi </math> and <math display="inline">l_\eta </math> are taken as the maximum projection of the velocities <math display="inline">u_\xi </math> and <math display="inline">u_\eta </math> along the element sides (for triangles) and the element diagonals (for quadrilaterals), i.e.
605
606
{| class="formulaSCP" style="width: 100%; text-align: left;" 
607
|-
608
| 
609
{| style="text-align: left; margin:auto;" 
610
|-
611
| style="text-align: center;" | <math>l_i =\max ({d}_j^T {u}_i)\quad ,\quad i=\xi ,\eta </math>
612
|}
613
| style="width: 5px;text-align: right;" |  (39a)
614
|}
615
616
with
617
618
{| class="formulaSCP" style="width: 100%; text-align: left;" 
619
|-
620
| 
621
{| style="text-align: left; margin:auto;" 
622
|-
623
| style="text-align: center;" | <math>\begin{array}{ll}j=1,2,3 \hbox{ (for triangles) and }\\ j=1,2 \hbox{ (for quadrilaterals)}\end{array}</math>
624
|}
625
| style="width: 5px;text-align: right;" |  (39b)
626
|}
627
628
In Eq.(39a) <math display="inline">{u}_\xi </math> and <math display="inline">{u}_\eta </math> contain the global components of the  velocity vectors <math display="inline">\vec u_\xi </math> and <math display="inline">\vec u_\eta </math>, respectively. For triangles <math display="inline">{d}_j</math> are the element side vectors, whereas for quadrilaterals <math display="inline">{d}_j</math> are the element diagonal vectors <span id='citeF-27'></span>[[#cite-27|27]].
629
630
The next step is to transform Eq.(34) to global axes <math display="inline">x,y</math>. The resulting equation is
631
632
{| class="formulaSCP" style="width: 100%; text-align: left;" 
633
|-
634
| 
635
{| style="text-align: left; margin:auto;" 
636
|-
637
| style="text-align: center;" | <math>-{u}^T {\boldsymbol \nabla }\phi +{\boldsymbol \nabla }^T {D}_G {\boldsymbol \nabla }\phi -s\phi=0 </math>
638
|}
639
|}
640
641
where the global diffusion matrix <math display="inline">{D}_G</math> is given by
642
643
{| class="formulaSCP" style="width: 100%; text-align: left;" 
644
|-
645
| 
646
{| style="text-align: left; margin:auto;" 
647
|-
648
| style="text-align: center;" | <math>{D}_G={D}+\bar{D}</math>
649
|}
650
| style="width: 5px;text-align: right;" | (41a)
651
|}
652
653
where the global balancing diffusion matrix <math display="inline">\bar{D}</math> is
654
655
{| class="formulaSCP" style="width: 100%; text-align: left;" 
656
|-
657
| 
658
{| style="text-align: left; margin:auto;" 
659
|-
660
| style="text-align: center;" | <math>\bar {D} = {T}^T \bar {D}^\prime {T}</math>
661
|}
662
| style="width: 5px;text-align: right;" | (41b)
663
|}
664
665
===Remark===
666
667
Similarly as for the 1D problems, a negative value of the parameters <math display="inline">\beta _\xi </math> and <math display="inline">\beta _\eta </math> indicates that no stabilization is needed along the directions <math display="inline">\xi </math> and <math display="inline">\eta </math>, respectively. A zero value of the corresponding stabilization parameter is chosen in this case.
668
669
===Remark===
670
671
The expressions of <math display="inline">\beta _\xi </math> and <math display="inline">\beta _\eta </math> in Eq.(37) can be simplified to
672
673
{| class="formulaSCP" style="width: 100%; text-align: left;" 
674
|-
675
| 
676
{| style="text-align: left; margin:auto;" 
677
|-
678
| style="text-align: center;" | <math>\begin{array}{l}\beta _\xi \simeq \beta _{\xi _c} =\left[\displaystyle{w_\xi \over 6} +|\gamma _\xi | -1\right]\\ \beta _\eta \simeq \beta _{\eta _c} =\left[\displaystyle{w_\eta \over 6} +|\gamma _\eta | -1\right] \end{array} </math>
679
|}
680
| style="width: 5px;text-align: right;" | (42)
681
|}
682
683
This avoids the computation of the sign of the solution and of its first and second derivatives. The expressions of <math display="inline">\beta _{\xi _c}</math> and <math display="inline">\beta _{\eta _c}</math> in Eq.(42) are equivalent to that of the 1D critical stabilization parameter <math display="inline">\beta _c</math> of Eq.(23). The main difference is that the computation of the local directions <math display="inline">\xi </math> and <math display="inline">\eta </math> is still mandatory in the multidimensional case and, therefore, the nonlinearity of the process can not be avoided here.
684
685
===5.1 Computation of the principal curvature axes for linear elements===
686
687
Excellent results have been obtained in our work ''by approximating the main  curvature direction <math>\vec \xi </math> by the direction of the gradient vector'' <math display="inline">{\boldsymbol \nabla } \phi </math>.
688
689
This simplification allows us to estimate the direction <math display="inline">\vec {\xi }</math> in a very economical manner as the gradient vector <math display="inline">{\boldsymbol \nabla } \bar \phi </math> can be directly computed at any point of a linear element. Direction <math display="inline">\vec {\eta }</math> is taken orthogonal to that of <math display="inline">\vec {\xi }</math> in an anti-clockwise sense.
690
691
For linear triangles <math display="inline">{\boldsymbol \nabla }\bar \phi </math> is constant within the element. For four node quadrilaterals <math display="inline">{\boldsymbol \nabla } \bar \phi </math> varies linearly. We have assumed in this case that the direction of <math display="inline">\vec \xi </math> is constant within the element and equal to the direction of vector <math display="inline">{\boldsymbol \nabla } \bar \phi </math> computed at the element center.
692
693
The computation of the signs of the second derivative of <math display="inline">\bar \phi </math> in Eq.(38) involves the following steps: 1) recovery of the cartesian first derivative field at the nodes using a nodal averaging procedure; 2) computation of the second derivative tensor at the element center and 3) transformation of this tensor to the local axes <math display="inline">\xi </math> and <math display="inline">\eta </math>.
694
695
We note that in problems where positive values of <math display="inline">\bar \phi </math> are prescribed at the Dirichlet boundary, the signs of <math display="inline">S_{\xi _2}</math>, <math display="inline">S_{\eta _2}</math> are positive due to the convexity of the numerical solution.
696
697
As mentioned above the dependence of the balancing diffusion matrix <math display="inline">\bar{D}</math> with the principal  curvature directions <math display="inline">\vec \xi </math> and <math display="inline">\vec {\eta }</math> introduces a nonlinearity in the solution process. A simple and effective iterative algorithm is described next.
698
699
===5.2 General iterative scheme===
700
701
A stabilized numerical solution can be found by the following algorithm.
702
703
'''Step 1''' . For elements in the interior of the domain choose <math display="inline">{}^1{\boldsymbol \xi } ={u}</math>, i.e. the gradient direction is taken coincident with the velocity direction. If <math display="inline">{u}=0</math> then <math display="inline">{}^1{\boldsymbol \xi }</math> is taken coincident with the global <math display="inline">{x}</math> axis.
704
705
In elements adjacent to a boundary choose <math display="inline">{}^1{\boldsymbol \xi } ={n }</math> where '''n ''' is the normal to the boundary.
706
707
Compute <math display="inline">{}^1{\boldsymbol \eta }, {}^1\bar  {D}'</math>, <math display="inline">{}^1\bar  {D}</math> and <math display="inline">{}^1{D}_G</math> using the expressions of <math display="inline">\beta _\xi </math> and <math display="inline">\beta _\eta </math> from Eq.(42).
708
709
Solve for <math display="inline">{}^1\bar{\boldsymbol \phi }</math>.
710
711
Verify that the solution is stable. This implies that there are not undershoots or overshoots in the numerical results with respect to the expected physical values. If the numerical solution is unstable, then go to step 2.
712
713
'''Step 2''' . For all elements, compute at the element center <math display="inline">{}^2{\boldsymbol \xi } = {\boldsymbol \nabla }^1\bar \phi </math>. Then compute <math display="inline">{}^2{\boldsymbol \eta } , {}^2\bar{D}'</math>, <math display="inline">{}^2\bar  {D}</math> and <math display="inline">{}^2{D}_G</math> using the expressions of <math display="inline">\beta _\xi </math> and <math display="inline">\beta _\eta </math> from Eqs.(37).
714
715
Solve for <math display="inline">{}^2\bar{\boldsymbol \phi }</math>.
716
717
'''Step 3''' . Estimate the convergence of the process. We have chosen the following convergence norm
718
719
{| class="formulaSCP" style="width: 100%; text-align: left;" 
720
|-
721
| 
722
{| style="text-align: left; margin:auto;" 
723
|-
724
| style="text-align: center;" | <math>\Vert \phi \Vert = {1\over N\bar \phi _{max}} \left[\sum \limits _{j=1}^n \left({}^{i}\bar \phi _j - {}^{i-1}\bar \phi _j \right)^2 \right]^{1/2} \le \varepsilon  </math>
725
|}
726
| style="width: 5px;text-align: right;" | (43)
727
|}
728
729
where <math display="inline">N</math> is the total number of nodes in the mesh and <math display="inline">\phi _{max}</math> is the maximum prescribed value at the Dirichlet boundary (if <math display="inline">\bar \phi _{max}=0</math> then <math display="inline">\bar  \phi _{max}=1</math>). In above steps the left upper indices denote the iteration number.
730
731
In the examples shown in the next section <math display="inline">\varepsilon =10^{-3}</math> has been taken.
732
733
If condition (43) is not satisfied, repeat steps 2 and 3 until convergence.
734
735
===Remark===
736
737
For the advective-diffusive problems (i.e. <math display="inline">s=0</math>) the  expression of the balancing diffusion matrix in the interior elements for step 1 coincides with the standard (linear) SUPG form <span id='citeF-27'></span>[[#cite-27|27]].
738
739
===Remark===
740
741
An  alternative solution scheme is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix.
742
743
==6 1D NUMERICAL EXAMPLES==
744
745
The examples presented in this section are solved in a 1D domain discretized with eight two-node linear elements of unit size. The examples are solved with the same Dirichlet conditions <math display="inline">\phi _1^p =8</math> and <math display="inline">\phi _9^p=3</math> and two different values of <math display="inline">\gamma </math> and <math display="inline">w</math> (<math display="inline">\gamma =1, w=20</math> and <math display="inline">\gamma =10, w=20</math>). We note that for both cases the instability condition (<math display="inline">\beta _c>0</math>) is violated and, hence, the Galerkin solution should yield non-physical results.
746
747
Figures 2 and 3 show the numerical results obtained with the standard  Galerkin method (<math display="inline">\beta =0</math>) and using the element (<math display="inline">\beta ^e</math>) and critical (<math display="inline">\beta _c</math>) values of the stabilization parameter  given by Eqs.(22) and (23), respectively. The exact analytical solution is also shown for comparison.
748
749
Table 1 shows the nodal values of the results of the example of Figure 3 for comparison with the 2D solutions presented in the next section.
750
751
The following conclusions are drawn from the 1D results.
752
753
<ol>
754
755
<li>The Galerkin solution (<math display="inline">\beta =0</math>) is unstable for both problems, as expected. </li>
756
757
<li>The solution obtained with the critical value <math display="inline">\beta _c</math> is always stable, although it yields slightly overdiffusive results in both cases. </li>
758
759
<li>The results obtained with <math display="inline">\beta ^e</math>  are less diffusive and more accurate than those obtained with <math display="inline">\beta _c</math>. The explanation is that the sign of the ratio <math display="inline">{S_1/S_2}</math> is negative in the region close to the left end point of the domain. This naturally reduces the value of the stabilizing diffusion parameter <math display="inline">\beta </math> in Eq.(21) with respect to that of <math display="inline">\beta _c</math> in Eq.(23) where the sign effect is not relevant.  </li>
760
761
<li>The FIC algorithm provides a stabilized solution for Dirichlet problems when there is a negative streamwise gradient of the solution. This is an advantage versus SUPG, GLS and SGS methods using a single stabilization parameter which fail in some cases for these type of problems <span id='citeF-12'></span>[[#cite-12|12]]. </li>
762
763
</ol>
764
765
Above conclusions have been confirmed in the solution of a wider collection of 1D problems presented in [18].
766
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
767
|-
768
|[[Image:draft_Samper_447243531-Fig2.png|600px|ϕ₁<sup>p</sup>=8, ϕ₉<sup>p</sup>=3, γ=1 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution.]]
769
|- style="text-align: center; font-size: 75%;"
770
| colspan="1" | <math>\phi _1^p=8, \phi _9^p =3, \gamma =1</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution.
771
|}
772
773
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
774
|-
775
|[[Image:draft_Samper_447243531-Fig6.png|600px|ϕ₁<sup>p</sup>=8, ϕ₉<sup>p</sup>=3, γ=10 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution.]]
776
|- style="text-align: center; font-size: 75%;"
777
| colspan="1" | <math>\phi _1^p =8, \phi _9^p =3, \gamma =10</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution.
778
|}
779
780
==7 2D EXAMPLES==
781
782
The analysis domain in the  first two 2D examples presented is a square of size 8 units. The problems are solved first with relatively coarse meshes of <math display="inline">8\times 8</math> four node bi-linear square elements and <math display="inline">8\times 8\times 2</math> linear triangles.
783
784
The boundary conditions for both examples are <math display="inline">\phi ^p =8</math> and <math display="inline">\phi ^p =3</math> at the boundaries <math display="inline">x=0</math> and <math display="inline">x=8</math>, respectively and zero normal flux at <math display="inline">y=0</math> and <math display="inline">y=8</math>. This reproduces the condition of the two 1D examples solved in the previous section. The first example is analized for <math display="inline">{u} = [2,0]^T</math>, <math display="inline">k=1</math> and <math display="inline">s=20</math> giving <math display="inline">w=20</math>, <math display="inline">\gamma _x=1</math> and <math display="inline">\gamma _y=0</math> which corresponds to the first 1D example (Figure 2). The correct solution for this problem has a boundary layer in the vecinity of the  two sides at <math display="inline">x=0</math> and <math display="inline">x=8</math> where <math display="inline">\phi </math> is prescribed (Figure 4). The numerical results obtained with the standard  Galerkin solution  are oscillatory as expected. The stabilized FIC formulation elliminates the oscillations and yields the correct physical solution. Good results are obtained for both meshes of linear rectangles and triangles (Figures 4 and 5).
785
786
Results labelled as FIC-1 and FIC-2 in the figures correspond to those obtained in the first and second iteration of the algorithm presented in Section 5.2, respectively. We note that the FIC-1 results agree precisely with those obtained in the 1D case for <math display="inline">\beta =\beta _c</math>, whereas the FIC-2 results agree with the more accurate 1D values obtained with the element stabilization parameter <math display="inline">\beta ^e</math> (see Figure 2).
787
788
The second example is similar to the first one with <math display="inline">{u} =[20,0]^T</math>, <math display="inline">k=1</math> and <math display="inline">s=20</math> giving <math display="inline">w=20</math>, <math display="inline">\gamma _x =10</math> and <math display="inline">\gamma _y =0</math>. These values correspond to the second 1D problem of the previous section (Figure 3). The Galerkin solution is again oscillatory, whereas the FIC results are physically sound (Figures 6 and 7). Once more  the FIC-1 and FIC-2 results  are in good agreement with the 1D values shown in Figure 3 for <math display="inline">\beta _c</math> and <math display="inline">\beta _e</math>, respectively for both meshes of square and triangular elements. The coincidence of the 1D and 2D results for this problem can be clearly seen in Table 1.
789
790
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
791
|-
792
|[[Image:draft_Samper_447243531-Figura4a.png|600px|]]
793
|-
794
|[[Image:draft_Samper_447243531-Figura4b.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [2,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=1 and γ<sub>y</sub>= 0. Galerkin and FIC solutions for a mesh of 8 ×8 four node square elements.]]
795
|- style="text-align: center; font-size: 75%;"
796
| colspan="1" | 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=8</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=8</math>. <math>{u} = [2,0]^T</math>, <math>k=1</math>, <math>s=20</math>, <math>w=20</math>, <math>\gamma _x=1</math> and <math>\gamma _y= 0</math>. Galerkin and FIC solutions for a mesh of <math>8 \times 8</math> four node square elements.
797
|}
798
799
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
800
|-
801
|[[Image:draft_Samper_447243531-Figura5a.png|600px|]]
802
|-
803
|[[Image:draft_Samper_447243531-Figura5b.png|600px|Solution of problem of Figure 4 with a mesh of 8 ×8×2 linear triangles.]]
804
|- style="text-align: center; font-size: 75%;"
805
| colspan="1" | Solution of problem of Figure 4 with a mesh of <math>8 \times 8\times 2</math> linear triangles.
806
|}
807
808
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
809
|-
810
|[[Image:draft_Samper_447243531-Figura6a.png|600px|]]
811
|-
812
|[[Image:draft_Samper_447243531-Figura6b.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [20,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=10 and γ<sub>y</sub>= 0. Galerkin and FIC solutions for a mesh of 8 ×8 four node square elements.]]
813
|- style="text-align: center; font-size: 75%;"
814
| colspan="1" | 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=8</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=8</math>. <math>{u} = [20,0]^T</math>, <math>k=1</math>, <math>s=20</math>, <math>w=20</math>, <math>\gamma _x=10</math> and <math>\gamma _y= 0</math>. Galerkin and FIC solutions for a mesh of <math>8 \times 8</math> four node square elements.
815
|}
816
817
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
818
|-
819
|[[Image:draft_Samper_447243531-Figura7a.png|600px|]]
820
|-
821
|[[Image:draft_Samper_447243531-Figura7b.png|600px|Solution of problem of Figure 5 with a mesh of 8 ×8×2 linear triangles.]]
822
|- style="text-align: center; font-size: 75%;"
823
| colspan="1" | Solution of problem of Figure 5 with a mesh of <math>8 \times 8\times 2</math> linear triangles.
824
|}
825
826
<br/>
827
828
829
{| class="wikitable" style="text-align: center; margin: 1em auto;"
830
|+ Table. 1 Comparison of 1D and 2D solutions for the advection-diffusion-absorption problem of Figure 3  (<math>\gamma _x =10</math>, <math>w=20</math>)
831
|- style="border-top: 2px solid;"
832
| colspan='5' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | '''1D'''
833
| colspan='4' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | '''2D (nodes along line A-A')'''
834
|- style="border-top: 2px solid;"
835
| colspan='5' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Figure 3
836
| colspan='2' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 4 node quads. (Fig. 6)
837
| colspan='2' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 3 node triangles (Fig. 7)
838
|- style="border-top: 2px solid;"
839
| style="border-left: 2px solid;border-right: 2px solid;" |  Node 
840
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\bar \phi (\beta =0</math>)
841
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\bar \phi (\beta ^e</math>)
842
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\bar \phi (\beta _c</math>)
843
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\phi </math>(exact)
844
| style="border-left: 2px solid;border-right: 2px solid;" | FIC-1
845
| style="border-left: 2px solid;border-right: 2px solid;" | FIC-2
846
| style="border-left: 2px solid;border-right: 2px solid;" | FIC-1
847
| style="border-left: 2px solid;border-right: 2px solid;" | FIC-2
848
|- style="border-top: 2px solid;"
849
| style="border-left: 2px solid;border-right: 2px solid;" |  1 
850
| style="border-left: 2px solid;border-right: 2px solid;" | 8,00 
851
| style="border-left: 2px solid;border-right: 2px solid;" | 8 
852
| style="border-left: 2px solid;border-right: 2px solid;" | 8 
853
| style="border-left: 2px solid;border-right: 2px solid;" | 8
854
| style="border-left: 2px solid;border-right: 2px solid;" | 8
855
| style="border-left: 2px solid;border-right: 2px solid;" | 8
856
| style="border-left: 2px solid;border-right: 2px solid;" | 8
857
| style="border-left: 2px solid;border-right: 2px solid;" | 8
858
|- style="border-top: 2px solid;"
859
| style="border-left: 2px solid;border-right: 2px solid;" |  2 
860
| style="border-left: 2px solid;border-right: 2px solid;" | 2,94 
861
| style="border-left: 2px solid;border-right: 2px solid;" | 3,06 
862
| style="border-left: 2px solid;border-right: 2px solid;" | 4 
863
| style="border-left: 2px solid;border-right: 2px solid;" | 3,08 
864
| style="border-left: 2px solid;border-right: 2px solid;" | 3,99 
865
| style="border-left: 2px solid;border-right: 2px solid;" | 3,057  
866
| style="border-left: 2px solid;border-right: 2px solid;" | 4,0 
867
| style="border-left: 2px solid;border-right: 2px solid;" | 3,059 
868
|- style="border-top: 2px solid;"
869
| style="border-left: 2px solid;border-right: 2px solid;" |  3 
870
| style="border-left: 2px solid;border-right: 2px solid;" | 1,32 
871
| style="border-left: 2px solid;border-right: 2px solid;" | 1,17 
872
| style="border-left: 2px solid;border-right: 2px solid;" | 2 
873
| style="border-left: 2px solid;border-right: 2px solid;" | 1,19 
874
| style="border-left: 2px solid;border-right: 2px solid;" | 2,00
875
| style="border-left: 2px solid;border-right: 2px solid;" | 1,170 
876
| style="border-left: 2px solid;border-right: 2px solid;" | 2,0 
877
| style="border-left: 2px solid;border-right: 2px solid;" | 1,167
878
|- style="border-top: 2px solid;"
879
| style="border-left: 2px solid;border-right: 2px solid;" |  4 
880
| style="border-left: 2px solid;border-right: 2px solid;" | 1,80 
881
| style="border-left: 2px solid;border-right: 2px solid;" | 0,447 
882
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
883
| style="border-left: 2px solid;border-right: 2px solid;" | 0,457 
884
| style="border-left: 2px solid;border-right: 2px solid;" | 1,00
885
| style="border-left: 2px solid;border-right: 2px solid;" | 0,448 
886
| style="border-left: 2px solid;border-right: 2px solid;" | 1,0 
887
| style="border-left: 2px solid;border-right: 2px solid;" | 0,452
888
|- style="border-top: 2px solid;"
889
| style="border-left: 2px solid;border-right: 2px solid;" |  5 
890
| style="border-left: 2px solid;border-right: 2px solid;" | 0,599 
891
| style="border-left: 2px solid;border-right: 2px solid;" | 0,172 
892
| style="border-left: 2px solid;border-right: 2px solid;" | 0,5 
893
| style="border-left: 2px solid;border-right: 2px solid;" | 0,176 
894
| style="border-left: 2px solid;border-right: 2px solid;" | 0,49 
895
| style="border-left: 2px solid;border-right: 2px solid;" | 0,172
896
| style="border-left: 2px solid;border-right: 2px solid;" | 0,499 
897
| style="border-left: 2px solid;border-right: 2px solid;" | 0,166
898
|- style="border-top: 2px solid;"
899
| style="border-left: 2px solid;border-right: 2px solid;" |  6 
900
| style="border-left: 2px solid;border-right: 2px solid;" | -0,633 
901
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0646 
902
| style="border-left: 2px solid;border-right: 2px solid;" | 0,25 
903
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0677 
904
| style="border-left: 2px solid;border-right: 2px solid;" | 0,248
905
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0648 
906
| style="border-left: 2px solid;border-right: 2px solid;" | 0,2501
907
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0681
908
|- style="border-top: 2px solid;"
909
| style="border-left: 2px solid;border-right: 2px solid;" |  7 
910
| style="border-left: 2px solid;border-right: 2px solid;" | 1,16 
911
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0264 
912
| style="border-left: 2px solid;border-right: 2px solid;" | 0,125 
913
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0261 
914
| style="border-left: 2px solid;border-right: 2px solid;" | 0,125
915
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0255 
916
| style="border-left: 2px solid;border-right: 2px solid;" | 0,1250
917
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0257
918
|- style="border-top: 2px solid;"
919
| style="border-left: 2px solid;border-right: 2px solid;" |  8 
920
| style="border-left: 2px solid;border-right: 2px solid;" | -1,83 
921
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0073
922
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0625
923
| style="border-left: 2px solid;border-right: 2px solid;" | 0,01
924
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0615 
925
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0101 
926
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0624
927
| style="border-left: 2px solid;border-right: 2px solid;" | 0,0072
928
|- style="border-top: 2px solid;border-bottom: 2px solid;"
929
| style="border-left: 2px solid;border-right: 2px solid;" |  9 
930
| style="border-left: 2px solid;border-right: 2px solid;" | 3
931
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
932
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
933
| style="border-left: 2px solid;border-right: 2px solid;" | 3
934
| style="border-left: 2px solid;border-right: 2px solid;" | 3
935
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
936
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
937
| style="border-left: 2px solid;border-right: 2px solid;" | 3
938
939
|}
940
941
Note that, similarly to the 1D case, the FIC-2 results are more accurate (less diffusive) than those obtained in the first iteration (FIC-1). This is due to the more precise evaluation of <math display="inline">\beta _s</math> and <math display="inline">\beta _\eta </math> in Eqs.(37) accounting for the correct sign of all the terms.
942
943
Figures 8&#8211;11 show results for the two 2D problems above described solved now with relatively coarse unstructured meshes of linear triangles and quadrilaterals. The effectiveness and accuracy of the FIC iterative scheme is again noticeable in all cases. Note the agreement of the FIC-2 results of Figures 10 and 11 with the exact solution for the equivalent 1D problem of Figure 3.
944
945
Figure 12 presents the solution of a similar problem where the values of <math display="inline">\phi </math> are prescribed at the four boundaries. The solution domain has now 10 units and the problem is solved first with a mesh of <math display="inline">10\times 10</math> four node square elements. Details of the physical parameters are given in Figure 12. Excellent results are again obtained with the FIC scheme. Similar good results are obtained with a structured mesh of linear triangles (Figure 13) as well as with  non structured meshes of linear quadrilateral and triangles (Figures 14 and 15).
946
947
The effectiveness of the FIC scheme for a diffusive-absorptive problem with Dirichlet boundary conditions is shown in Figure 16. The results shown have been obtained with structured meshes of linear quadrilateral and triangles. Note that the four boundary layers  are well captured in the first step of the iterative solution. Similar good results have also been obtained with unstructured meshes not shown here.
948
949
The final example is a standard benchmark problem of advection-diffusion where sharp layers appear at both the boundary and the interior of the domain. The problem is the advective-diffusive transport of <math display="inline">\phi </math> in a square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source terms (i.e. <math display="inline">Q=0</math> and <math display="inline">s=0</math>). Figure 17 displays the SUPG solution and FIC results obtained after two iterations using a structured mesh of <math display="inline">20\times 20</math> linear four node square elements. It is remarkable that the FIC results capture the sharp gradient zones at the boundaries where <math display="inline">\phi </math> is prescribed to zero and at the interior of the domain and elliminate all the spurious oscillations present in the SUPG method.
950
951
Similar good results obtained with the FIC method for a wide range of advective-diffusive problems are presented in <span id='citeF-27'></span>[[#cite-27|27]]. Recent applications of the FIC method to incompressible fluid flow problems are reported in [37].
952
953
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
954
|-
955
|[[Image:draft_Samper_447243531-Figura8a.png|600px|]]
956
|-
957
|[[Image:draft_Samper_447243531-Figura8b.png|600px|Solution of problem of Figure 4 with an unstructured mesh of 209 four node bi-linear quadrilaterals]]
958
|- style="text-align: center; font-size: 75%;"
959
| colspan="1" | Solution of problem of Figure 4 with an unstructured mesh of 209 four node bi-linear quadrilaterals
960
|}
961
962
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
963
|-
964
|[[Image:draft_Samper_447243531-Figura9a.png|600px|]]
965
|-
966
|[[Image:draft_Samper_447243531-Figura9b.png|600px|Solution of problem of Figure 4 with an unstructured mesh of 176 three node linear triangles]]
967
|- style="text-align: center; font-size: 75%;"
968
| colspan="1" | Solution of problem of Figure 4 with an unstructured mesh of 176 three node linear triangles
969
|}
970
971
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
972
|-
973
|[[Image:draft_Samper_447243531-Figura10a.png|600px|]]
974
|-
975
|[[Image:draft_Samper_447243531-Figura10b.png|600px|Solution of problem of Figure 6 with an unstructured mesh of 209 four node bi-linear quadrilaterals]]
976
|- style="text-align: center; font-size: 75%;"
977
| colspan="1" | Solution of problem of Figure 6 with an unstructured mesh of 209 four node bi-linear quadrilaterals
978
|}
979
980
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
981
|-
982
|[[Image:draft_Samper_447243531-Figura11a.png|600px|]]
983
|-
984
|[[Image:draft_Samper_447243531-Figura11b.png|600px|Solution of problem of Figure 6 with an unstructured mesh of 176 three node triangles]]
985
|- style="text-align: center; font-size: 75%;"
986
| colspan="1" | Solution of problem of Figure 6 with an unstructured mesh of 176 three node triangles
987
|}
988
989
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
990
|-
991
|[[Image:draft_Samper_447243531-Figura12a.png|600px|]]
992
|-
993
|[[Image:draft_Samper_447243531-Figura12b.png|600px|2D advection-diffusion-absorption problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=50 along x=0 and y=10 and ϕ<sup>p</sup>=100 along x=10 and y=0, u= [0,15]<sup>T</sup>, k=1, s=12, w=12, γₓ=0 and γ<sub>y</sub>= 7.5. Galerkin and FIC solutions for a mesh of 10 ×10 four node bi-linear square elements.]]
994
|- style="text-align: center; font-size: 75%;"
995
| colspan="1" | 2D advection-diffusion-absorption problem over a square domain of size equal to 10 units. <math>\phi ^p=50</math> along <math>x=0</math> and <math>y=10</math> and <math>\phi ^p=100</math> along <math>x=10</math> and <math>y=0</math>, <math>{u}= [0,15]^T</math>, <math>k=1</math>, <math>s=12</math>, <math>w=12</math>, <math>\gamma _x=0</math> and <math>\gamma _y= 7.5</math>. Galerkin and FIC solutions for a mesh of <math>10 \times 10</math> four node bi-linear square elements.
996
|}
997
998
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
999
|-
1000
|[[Image:draft_Samper_447243531-Figura13a.png|600px|]]
1001
|-
1002
|[[Image:draft_Samper_447243531-Figura13b.png|600px|Solution of the problem of Figure 12 with an unstructured mesh of 432 four node bi-linear quadrilaterals]]
1003
|- style="text-align: center; font-size: 75%;"
1004
| colspan="1" | Solution of the problem of Figure 12 with an unstructured mesh of 432 four node bi-linear quadrilaterals
1005
|}
1006
1007
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1008
|-
1009
|[[Image:draft_Samper_447243531-Figura14a.png|600px|]]
1010
|-
1011
|[[Image:draft_Samper_447243531-Figura14b.png|600px|Solution of problem of Figure 12 with an structured mesh of 10×10×2  three node linear triangles]]
1012
|- style="text-align: center; font-size: 75%;"
1013
| colspan="1" | Solution of problem of Figure 12 with an structured mesh of <math>10\times 10\times 2</math>  three node linear triangles
1014
|}
1015
1016
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1017
|-
1018
|[[Image:draft_Samper_447243531-Figura15a.png|600px|]]
1019
|-
1020
|[[Image:draft_Samper_447243531-Figura15b.png|600px|Solution of problem of Figure 12 with an unstructured mesh of 780 three node triangles]]
1021
|- style="text-align: center; font-size: 75%;"
1022
| colspan="1" | Solution of problem of Figure 12 with an unstructured mesh of 780 three node triangles
1023
|}
1024
1025
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1026
|-
1027
|[[Image:draft_Samper_447243531-Figura16a.png|600px|]]
1028
|-
1029
|[[Image:draft_Samper_447243531-Figura16b.png|600px|Diffusive-absorptive problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=50 over x=0 and y=10 and ϕ<sup>p</sup>=100 over x=10 and y=0, u= [0,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=0 and γ<sub>y</sub>= 0. Galerkin and FIC solutions obtained with structured  meshes of four node quadrilaterals and linear triangles.]]
1030
|- style="text-align: center; font-size: 75%;"
1031
| colspan="1" | Diffusive-absorptive problem over a square domain of size equal to 10 units. <math>\phi ^p=50</math> over <math>x=0</math> and <math>y=10</math> and <math>\phi ^p=100</math> over <math>x=10</math> and <math>y=0</math>, <math>{u}= [0,0]^T</math>, <math>k=1</math>, <math>s=20</math>, <math>w=20</math>, <math>\gamma _x=0</math> and <math>\gamma _y= 0</math>. Galerkin and FIC solutions obtained with structured  meshes of four node quadrilaterals and linear triangles.
1032
|}
1033
1034
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1035
|-
1036
|[[Image:draft_Samper_447243531-Figura16c.png|600px|]]
1037
|-
1038
|[[Image:draft_Samper_447243531-Figura16d.png|600px|(cont.)]]
1039
|- style="text-align: center; font-size: 75%;"
1040
| colspan="1" | (cont.)
1041
|}
1042
1043
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1044
|-
1045
|[[Image:draft_Samper_447243531-Figure7.png|600px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 20×20 linear four node square elements]]
1046
|- style="text-align: center; font-size: 75%;"
1047
| colspan="1" | Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>20\times 20</math> linear four node square elements
1048
|}
1049
1050
==8 CONCLUSIONS==
1051
1052
The FIC-FEM formulation presented allows to obtain a stabilized and accurate solution for the advection-diffusion-absorption equation. For the 1D problem the formulation is equivalent to adding a non-linear diffusion term to the standard discretized equations which is  governed by a single stabilization parameter. The use of the constant critical value of the 1D stabilization parameter  provides a stabilized numerical solution in ''a single step''. A more accurate (less diffusive) solution can be obtained using the two step iterative scheme proposed.
1053
1054
The equivalence of the FIC  method with a nonlinear stabilizing diffusion term extends naturally to multidimensional problems using structured and unstructured meshes. The key step is to express the governing equations of the FIC formulation in the principal curvature directions of the solution. The resulting FIC equation is equivalent to adding a nonlinear diffusion matrix to the infinitessimal governing equations. The solution process becomes non linear and a simple iterative algorithm has been presented. The approximation of the main principal curvature direction by that of the gradient vector simplifies the computations in the iterative scheme. Excellent results have been obtained for all the 2D problems solved in just two iterations for structured and nonstructured meshes.
1055
1056
It is remarkable that, similarly to the 1D case, good stabilized results are obtained in the first iteration of the scheme proposed (FIC-1 results) and this may be sufficient for many practical cases. More accurate (less diffusive) results are obtained by performing a  second iteration at a relatively small additional computational cost.
1057
1058
==ACKNOWLEDGEMENTS==
1059
1060
The authors also thank Profs. C. Felippa and S.R. Idelsohn for many useful discussions.
1061
1062
This work has been sponsored by the ''Ministerio de Educación y Ciencia of Spain''.  Plan Nacional, Project numbers: DPI2001-2193, BIA2003-09078-C02-01, and DPI2004-07410-C03-02.
1063
1064
===BIBLIOGRAPHY===
1065
1066
<div id="cite-1"></div>
1067
'''[[#citeF-1|[1]]]'''  Zienkiewicz O.C.  and  Taylor R.L. ''The Finite Element Method''. Volume 3. 5th Edition, Butterworth-Heinemann, 2001.
1068
1069
<div id="cite-2"></div>
1070
'''[[#citeF-2|[2]]]''' Tezduyar T.E. and Park Y.J. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 307&#8211;325, 1986.
1071
1072
<div id="cite-3"></div>
1073
'''[3]''' Tezduyar T.E., Park Y.J. and Deans H.A. Finite element procedures for time-dependent convection-diffusion-reaction systems. ''Int. J. Num. Meth. Fluids'', '''7''', 1013&#8211;1033, 1987.
1074
1075
<div id="cite-4"></div>
1076
'''[[#citeF-4|[4]]]'''  Idelsohn S., Nigro N. Storti M. and Buscaglia G. A Petrov-Galerkin formulation for advection-reaction-diffusion problems. ''Comput. Methods Appl. Mech. Engrg.'', '''136''', 27&#8211;46, 1996.
1077
1078
<div id="cite-5"></div>
1079
'''[[#citeF-5|[5]]]'''  Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', '''156''', 186&#8211;210, 1998.
1080
1081
<div id="cite-6"></div>
1082
'''[6]'''  Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. ''Comput. Meth. Appl. Mech. Engrg.'', '''188''', 61&#8211;82, 2000.
1083
1084
<div id="cite-7"></div>
1085
'''[7]'''  Burman E. and Ern A. Nonlinear diffusion and discrete maxium principle for stabilized Galerkin approximations of the convection-diffusion-reaction equation. ''Comput. Meth. Appl. Mech. Engrg.'', '''191''', 3833&#8211;3855, 2002.
1086
1087
<div id="cite-8"></div>
1088
'''[[#citeF-8|[8]]]'''  Harari I. and  Hughes T.J.R. Finite element methods for the Helmholtz equation in an exterior domain: model problems. ''Comput. Meth. Appl. Mech. Engrg.'', '''87''', 59&#8211;96, 1991.
1089
1090
<div id="cite-9"></div>
1091
'''[9]'''   Harari I. and  Hughes T.J.R. Galerkin/least squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. ''Comput. Meth. Appl. Mech. Engrg.'', '''98''', 411&#8211;454, 1992.
1092
1093
<div id="cite-10"></div>
1094
'''[10]'''   Harari I. and  Hughes T.J.R. Stabilized finite element method for steady advection-diffusion with production. ''Comput. Meth. Appl. Mech. Engrg.'', '''115''', 165&#8211;191, 1994.
1095
1096
<div id="cite-11"></div>
1097
'''[11]'''  Harari I., Grosh K., Hughes T.J.R., Malhotra M.,  Pinsky P.M., Stewart J.R. and  Thompson L.L. Recent development in finite element methods for structural acoustics. ''Archives of Computational Mechanics in Engineering'', '''3, 2-3''', 131&#8211;309, 1996.
1098
1099
<div id="cite-12"></div>
1100
'''[[#citeF-12|[12]]]'''  Hauke G. and Garcia-Olivares A. Variational subgrid scale formulations for the advection-diffusion-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'' 2000; '''190''':6847&#8211;6865.
1101
1102
<div id="cite-13"></div>
1103
'''[13]''' Hauke  G. A simple subgrid scale stabilized method for the advection-diffusion reaction equation. ''Comput. Methods Appl. Mech. Engrg.'' 2002; '''191''':2925&#8211;2947.
1104
1105
<div id="cite-14"></div>
1106
'''[[#citeF-14|[14]]]'''  Brezzi F., Hauke G., Marin L.D. and Sangalli S. Link-cutting bubbles for the stabilization of convection-diffusion-reaction problems. ''Mathematical Models and Methods in Applied Sciences'', World Scientific Publishing Company, 2002.
1107
1108
<div id="cite-15"></div>
1109
'''[[#citeF-15|[15]]]'''  Felippa C.A. and  Oñate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Submitted to ''Comput. Mechanics''.
1110
1111
<div id="cite-16"></div>
1112
'''[16]'''  Idelsohn S.R.,  Heinrich J.C. and Oñate E. Petrov-Galerkin methods for the transient advective-diffusive equation with sharp gradients. ''Int. J. Num. Meth. Engng.'', '''39''', 1455&#8211;1473, 1996.
1113
1114
<div id="cite-17"></div>
1115
'''[17]'''  Harari I. Stability of semidiscrete advection-diffusion in transient computation. Proceedings of ''6th World Congress on Computational Mechanics'', Beijing, Sept. 2004, Z.H. Yao, M.W. Yuan and W.X. Zhong (Eds.), Tsinghua Univ. Press-Springer.
1116
1117
<div id="cite-18"></div>
1118
'''[18]''' Oñate E., Miquel, J. and Hauke, G. Stabilized formulation for the advection-diffusion-reaction equations using finite calculus and linear finite elements. Submittted in ''Comput. Methods Appl. Mech. Engrg.'', March 2005.
1119
1120
<div id="cite-19"></div>
1121
'''[19]''' 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.
1122
1123
<div id="cite-20"></div>
1124
'''[20]'''  Oñate E. Possibilities of finite calculus in computational mechanics. ''Int. J. Num. Meth. Engng.'', '''60''', 255&#8211;281, 2004.
1125
1126
<div id="cite-21"></div>
1127
'''[21]''' 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.
1128
1129
<div id="cite-22"></div>
1130
'''[22]'''  Oñate E., Idelsohn S.R., Zienkiewicz O.C. and Taylor R.L. A finite point method in computational mechanics. Applications to convective transport and fluid flow. ''Int. J. Num. Meth. Engng.'', '''39''', 3839&#8211;3866, 1996.
1131
1132
<div id="cite-23"></div>
1133
'''[23]'''  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.
1134
1135
<div id="cite-24"></div>
1136
'''[24]'''  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.
1137
1138
<div id="cite-25"></div>
1139
'''[25]'''  Idelsohn S.R., Oñate E. and Del Pin F. A Lagrangian meshless finite element method applied to fluid-structure interaction problems. ''Computers and Structures'', '''81''', 655&#8211;671, 2003.
1140
1141
<div id="cite-26"></div>
1142
'''[26]'''   Oñate E., Idelsohn S.R., Del Pin F. and  Aubry R. The particle finite element method. An overview (PFEM). ''Int. J. Comput. Methods'', '''1''' ('''2'''), 267&#8211;307, 2004.
1143
1144
<div id="cite-27"></div>
1145
'''[[#citeF-27|[27]]]'''   Oñate E., García J. and Idelsohn S.R. Ship Hydrodynamics. ''Encyclopedia of Computational Mechanics'', T. Hughes, R. de Borst and E. Stein (Eds.), Vol. 3, Chapter 18, 579&#8211;607, J. Wiley, 2004.
1146
1147
<div id="cite-28"></div>
1148
'''[28]'''  Oñate E., Taylor R.L., Zienkiewicz O.C. and  Rojek J. A residual correction method based on finite calculus. ''Engineering Computations'', '''20''' (5/6), 629&#8211;658, 2003.
1149
1150
<div id="cite-29"></div>
1151
'''[29]''' Oñate E., Rojek J., Taylor R.L. and Zienkiewicz O.C. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. ''Int. J. Num. Meth. Engng.'', '''59''' (11), 1473&#8211;1500, 2004.
1152
1153
<div id="cite-30"></div>
1154
'''[30]''' Oñate E., Zárate F. and Idelsohn S.R. Finite element formulation for convection-diffusion problems with sharp gradients using finite calculus. ''Comput. Meth. Appl. Mech. Engng.'' Submitted Nov. 2004.
1155
1156
<div id="cite-31"></div>
1157
'''[31]'''  Hughes T.J.R. and 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.
1158
1159
<div id="cite-32"></div>
1160
'''[32]'''  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.
1161
1162
<div id="cite-33"></div>
1163
'''[33]'''  Hughes T.J.R.,  Mallet M. and Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. ''Comput. Methods Appl. Mech. Engrg.'', '''54''', 341&#8211;355, 1986.
1164
1165
<div id="cite-34"></div>
1166
'''[34]'''  Galeo A.C. and  Dutra do Carmo E.G. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. ''Comput. Methods Appl. Mech. Engrg.'', '''68''', 83&#8211;95, 1988.
1167
1168
<div id="cite-35"></div>
1169
'''[35]'''  Tezduyar T.E. Adaptive determination of the finite element stabilization parameters. ''Proceedings of the ECCOMAS Computational Fluid Dynamics Conference 2001'', Swansea,  Wales, UK, CD-Rom, 2001.
1170
1171
<div id="cite-36"></div>
1172
'''[36]'''  Tezduyar T.E. Computation of moving boundaries and interfaces and stabilization parameters, ''International Journal for Numerical Methods in Fluid'', '''43''', 555-575, 2003.
1173
1174
<div id="cite-37"></div>
1175
'''[37]''' Oñate E., García J., Idelsohn S. and Del Pin F. FIC formulation for finite element analysis of incompressible flows. Eulerian, Lagrangian and ALE approaches. Accepted for publication in ''Computational Methods in Applied Mechanics and Engineering'', 2005.
1176

Return to Onate et al 2004a.

Back to Top