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
2
3
==Abstract==
4
5
The paper describes some recent developments  in finite element and particle methods for analysis  of a wide range of bulk forming  processes. The developments include new stabilized linear triangles and tetrahedra  using finite calculus and a new procedure combining particle methods and finite element methods. Applications of the new numerical methods  to casting, forging and  other bulk metal forming problems and mixing processes are shown.
6
7
'''Keywords''': Bulk forming processes, stabilized finite element method, particle method, particle finite element method, mixing processes.
8
9
==1 INTRODUCTION==
10
11
The development of efficient and robust numerical methods for analysis of bulk forming problems  has been a subject of intensive research in recent years <span id=‘citeF-1’></span>[[#cite-1|[1-7]]]. Many of these problems require the solution of incompressible fluid flow situations (such as in mould filling problems) whereas in other cases (such as forging, rolling, extrusion, etc.) the numerical method must be able to account for the quasi/fully incompressible behaviour induced by the large plastic deformation. The solution of these problems has motivated the development of the so called ''stabilized numerical methods'' overcoming the two main sources of instability in the analysis of incompressible  continua, namely those originated by the high values of the convective terms in fluid flow situtions and those induced by the difficulty in satisfying the incompressibility condition.
12
13
Different approaches to solve both type of  problems in the context of the finite element method (FEM) have been recently developed <span id=‘citeF-8’></span>[[#cite-8|[8]]]. Traditionally, the underdiffusive character of the Galerkin FEM for high convection flows  has been corrected by adding some kind of artificial viscosity terms to the standard Galerkin equations <span id=‘citeF-8’></span>[[#cite-8|[8,9]]].
14
15
A popular way to overcome the problems with the incompressibility constraint in the FEM is by introducing a pseudo-compressibility in the continuum  and using implicit and explicit algorithms ''ad hoc'' such as artificial compressibility schemes <span id=‘citeF-10’></span>[[#cite-10|[10]]] and preconditioning techniques <span id=‘citeF-11’></span>[[#cite-11|[11]]]. Other FEM schemes with good stabilization properties for the convective and incompressibility terms  in fluid flows are based in Petrov-Galerkin (PG) techniques. The background of PG methods are the non-centred (upwind) schemes for computing the first derivatives of the convective operator in FD and FV methods <span id=‘citeF-8’></span>[[#cite-8|[8,9,12]]]. A general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many FEM of this kind  we can name the Streamline Upwind Petrov Galerkin (SUPG) method <span id=‘citeF-8’></span>[[#cite-8|[8,13-16]]], the Galerkin Least Square (GLS) method <span id=‘citeF-8’></span>[[#cite-8|[8,17]]], the Taylor-Galerkin method <span id='citeF-18'></span>[[#cite-18|[18]]], the Characteristic Galerkin method <span id='citeF-19'></span>[[#cite-19|[19]]] and its variant the Characteristic Based Split (CBS) method <span id='citeF-20'></span>[[#cite-20|[20,21]]], the pressure gradient operator method <span id='citeF-22'></span>[[#cite-22|[22]]] and the Subgrid Scale (SS) method <span id='citeF-23'></span>[[#cite-23|[23]]]. A good review of these methods can be found in <span id='citeF-24'></span>[[#cite-24|[24]]]. Extensions of the CBS and SS methods to treat incompressible problems in solid mechanics are reported in <span id='citeF-25'></span>[[#cite-25|[25,26,58]]] and <span id='citeF-32'></span>[[#cite-32|[27-29]]], respectively.
16
17
In this paper a different class of stabilized FEM for quasi and fully incompressible fluid and solid materials applicable to a wide range of bulk forming problems is presented. The starting point is  the modified governing differential equations of the continuous  problem formulated via a finite calculus (FIC) approach <span id='citeF-35'></span>[[#cite-35|[30,31]]]. The FIC method is based in invoking the classical balance (or equilibrium) laws in a  domain of ''finite size''. This introduces naturally additional terms in the  differential equations of infinitesimal continuum  mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin FEM. One of the main advantages of the FIC formulation versus other alternative approaches (such as mixed FEM, etc.) is that it allows  to solve incompressible fluid problems using  low order finite elements (such as linear triangles and tetrahedra) with equal order approximations for the velocity and pressure variables <span id='citeF-32'></span>[[#cite-32|[32-35]]]. The FIC formulation has been successfully used for analysis of fully or quasi incompressible  solids <span id='citeF-36'></span>[[#cite-36|[36,37]]].
18
19
The layout of the paper is the following. In the next section the basic FIC equations for incompressible flow problems formulated in an Eulerian frame are presented. The finite element discretization is introduced and the resulting discretized equations are detailed. A fractional step scheme for the transient solution is presented.
20
21
The stabilized Eulerian formulation is extended to account for thermal effects and the transport of the free surface which are needed for mould filling processes.
22
23
The following sections outline the  FIC formulation for analysis of quasi/fully incompressible solids using a  Lagrangian description and linear triangles and tetrahedra. An explicit algorithm for integrating in time the equations of motion of elasto-plastic solids in large-strain problems involving frictional contact  is described. Examples of application to a casting problem and some bulk forming processes are presented.
24
25
In the last part of the paper  a  Lagrangian formulation for fluid flow analysis is presented  as a straightforward extension of the formulation for solid mechanics. The procedure, called the ''Particle Finite Element Method'' (PFEM) <span id='citeF-38'></span>[[#cite-38|[38-40]]], treats the mesh nodes in the fluid and solid domains as dimensionless particles which can freely move,  an even separate from the main fluid domain, representing, for instance, the effect of liquid drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The main  advantage of the Lagrangian flow formulation is that the convective terms do not enter in the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes.  The final examples show that the PFEM is  a promising method to solve mould filling and casting problems, material mixing processes and many other bulk metal forming problems involving the interaction between solids and fluids which can be treated with the same formulation.
26
27
==2 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOWS==
28
29
The FIC governing equations for a viscous incompressible fluid can be written in an Eulerian frame of reference as <span id='citeF-30'></span>[[#cite-30|[30-34]]]
30
31
'''Momentum'''
32
33
{| class="formulaSCP" style="width: 100%; text-align: left;" 
34
|-
35
| 
36
{| style="text-align: left; margin:auto;width: 100%;" 
37
|-
38
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}=0 \quad \hbox{in }\Omega  </math>
39
|}
40
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
41
|}
42
43
'''Mass balance'''
44
45
{| class="formulaSCP" style="width: 100%; text-align: left;" 
46
|-
47
| 
48
{| style="text-align: left; margin:auto;width: 100%;" 
49
|-
50
| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_j {\partial r_d \over \partial x_j}}{1\over 2} h_j {\partial r_d \over \partial x_j}=0 \quad \hbox{in }\Omega  </math>
51
|}
52
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
53
|}
54
55
where
56
57
{| class="formulaSCP" style="width: 100%; text-align: left;" 
58
|-
59
| 
60
{| style="text-align: left; margin:auto;width: 100%;" 
61
|-
62
| 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>
63
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
64
|-
65
| style="text-align: center;" | <math> r_d = {\partial v_i \over \partial x_i}\qquad i,j = 1, n_d </math>
66
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
67
|}
68
|}
69
70
Above <math display="inline">\Omega </math> is the analysis domain which can evolve with time, <math display="inline">n_d</math> is the number of space dimensions (<math display="inline">n_d=3</math> for 3D problems), <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 the body forces and <math display="inline">s_{ij}</math> are the  deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression
71
72
{| class="formulaSCP" style="width: 100%; text-align: left;" 
73
|-
74
| 
75
{| style="text-align: left; margin:auto;width: 100%;" 
76
|-
77
| 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>
78
|}
79
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
80
|}
81
82
where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are
83
84
{| class="formulaSCP" style="width: 100%; text-align: left;" 
85
|-
86
| 
87
{| style="text-align: left; margin:auto;width: 100%;" 
88
|-
89
| 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>
90
|}
91
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
92
|}
93
94
The  boundary conditions  are written in the FIC approach as
95
96
{| class="formulaSCP" style="width: 100%; text-align: left;" 
97
|-
98
| 
99
{| style="text-align: left; margin:auto;width: 100%;" 
100
|-
101
| style="text-align: center;" | <math>n_j \sigma _{ij} -t_i^p + \underline{{1\over 2} h_j n_j r_{m_i}}{1\over 2} h_j n_j r_{m_i}=0 \quad \hbox{on }\Gamma _t </math>
102
|}
103
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
104
|}
105
106
{| class="formulaSCP" style="width: 100%; text-align: left;" 
107
|-
108
| 
109
{| style="text-align: left; margin:auto;width: 100%;" 
110
|-
111
| style="text-align: center;" | <math>v_j - v_j^p =0 \quad \hbox{on }\Gamma _u </math>
112
|}
113
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
114
|}
115
116
and the initial condition is <math display="inline">v_j =v_j^0</math> for <math display="inline">t=t_0</math>.
117
118
Summation convention for repeated indexes in products and derivatives is used unless otherwise specified.
119
120
In Eqs.(7) and (8) <math display="inline">t_i^p</math> and <math display="inline">v_j^p</math> are surface tractions and prescribed velocities on the boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _u</math>, respectively, <math display="inline">n_j</math> are the components of the unit normal vector to the boundary and <math display="inline">\sigma _{ij}</math> are the total stresses given by <math display="inline">\sigma _{ij}=s_{ij}-\delta _{ij}p</math>.
121
122
Eqs. (1) and (2) are obtained by invoking the classical balance equations in fluid mechanics in a domain of finite size and retaining higher order terms <span id='citeF-30'></span>[[#cite-30|[30]]]. The <math display="inline">h_i's</math> in above equations are characteristic lengths of the domain where the balance of momentum and mass is enforced. In Eq.(7) these lengths define the domain where equilibrium of boundary tractions is established. In the discretized problem the <math display="inline">h_{i}</math> coincide with a typical element dimension, as described in Section 5. Note that by making <math display="inline">h_i =0</math> in these equations the standard infinitesimal form of the fluid mechanics equations is recovered <span id='citeF-8'></span>[[#cite-8|[8,9,24]]].
123
124
Eqs.(1)&#8211;(8) are the starting  point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations using equal order interpolation for the velocity and pressure variables <span id='citeF-32'></span>[[#cite-32|[32-35]]]. Application of the FIC formulation to meshless analysis of fluid flow problems using the finite point method can be found in <span id='citeF-42'></span>[[#cite-42|[42]]].
125
126
'''Remark'''. In most metal forming processes the viscosity <math display="inline">\mu </math> is a non linear function of the strain rate and the yield stress of the material (which also depends on the temperature) <span id='citeF-8'></span>[[#cite-8|[8,43,44]]]. This dependence adds another non linearity to the problem.
127
128
'''Remark'''. In Eqs.(1)&#8211;(7) <math display="inline">\Omega </math> and <math display="inline">\Gamma _t</math> respectively denote the actual volume and the boundary of the domain where the governing equations are solved at each instant of the forming processes. This domain is considered here to be fixed in space (Eulerian approach). Moving free surfaces, such as in the case of mold filling processes, can be modelled by using standard level set and volume of fluid (VOF) techniques <span id='citeF-47'></span>[[#cite-47|[47,48]]]. An alternative is to use an arbitrary Lagragian-Eulerian (ALE) description or even a fully Lagrangian description to follow the motion of the particles during the forming process.  The Lagrangian description is typical in the case of solid mechanics problems (Section 5). A particular Lagrangian formulation for fluid flow problems involving large motion of the free surface, named the Particle Finite Element Method, is presented in Section 8.
129
130
===2.1 Stabilized integral forms===
131
132
From the momentum equations it can be obtained <span id='citeF-32'></span>[[#cite-32|[32-34]]]
133
134
{| class="formulaSCP" style="width: 100%; text-align: left;" 
135
|-
136
| 
137
{| style="text-align: left; margin:auto;width: 100%;" 
138
|-
139
| style="text-align: center;" | <math>{\partial r_d \over \partial x_i}\simeq {h_j\over 2a_i} {\partial r_{m_i} \over \partial x_j}\quad ,\quad \hbox{no sum in }i </math>
140
|}
141
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
142
|}
143
144
where
145
146
{| class="formulaSCP" style="width: 100%; text-align: left;" 
147
|-
148
| 
149
{| style="text-align: left; margin:auto;width: 100%;" 
150
|-
151
| style="text-align: center;" | <math>a_i = {2\mu \over 3} +{\rho v_i h_i\over 2}\quad ,\quad \hbox{no sum in }i </math>
152
|}
153
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
154
|}
155
156
Substituting Eq.(9) into Eq.(2) and retaining the terms involving the derivatives of <math display="inline">r_{m_i}</math> with respect to <math display="inline">x_i</math> only, leads to the following expression for the stabilized mass balance equation
157
158
{| class="formulaSCP" style="width: 100%; text-align: left;" 
159
|-
160
| 
161
{| style="text-align: left; margin:auto;width: 100%;" 
162
|-
163
| style="text-align: center;" | <math>r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0 </math>
164
|}
165
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
166
|}
167
168
with
169
170
{| class="formulaSCP" style="width: 100%; text-align: left;" 
171
|-
172
| 
173
{| style="text-align: left; margin:auto;width: 100%;" 
174
|-
175
| style="text-align: center;" | <math>\tau _i = \left({8\mu \over 3h_i^2}+{2\rho v_i\over h_i}\right)^{-1} </math>
176
|}
177
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
178
|}
179
180
The <math display="inline">\tau _i</math>'s in Eq.(12), when scaled by the density, are termed in the stabilization literature ''intrinsic time parameters''.
181
182
The weighted residual form of the momentum and mass balance equations (Eqs.(1) and (11)) is written as
183
184
{| class="formulaSCP" style="width: 100%; text-align: left;" 
185
|-
186
| 
187
{| style="text-align: left; margin:auto;width: 100%;" 
188
|-
189
| style="text-align: center;" | <math>\int _\Omega \delta v_i \left[r_{m_i} - {h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]d\Omega + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i^p + {h_j \over 2} n_j r_{m_i}) d\Gamma =0 </math>
190
|}
191
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
192
|}
193
194
{| class="formulaSCP" style="width: 100%; text-align: left;" 
195
|-
196
| 
197
{| style="text-align: left; margin:auto;width: 100%;" 
198
|-
199
| style="text-align: center;" | <math>\int _\Omega q \left[r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]d\Omega =0 </math>
200
|}
201
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
202
|}
203
204
where <math display="inline">\delta v_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocities and virtual pressure fields. Integrating by parts the <math display="inline">r_{m_i}</math> terms  leads to
205
206
{| class="formulaSCP" style="width: 100%; text-align: left;" 
207
|-
208
| 
209
{| style="text-align: left; margin:auto;width: 100%;" 
210
|-
211
| style="text-align: center;" | <math>\int _\Omega \delta v_i r_{m_i}d\Omega  + \int _{\Gamma _t} \delta v_i (\sigma _{ij} n_j - t_i)d\Gamma + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
212
|}
213
| style="width: 5px;text-align: right;white-space: nowrap;" |  (15a)
214
|}
215
216
{| class="formulaSCP" style="width: 100%; text-align: left;" 
217
|-
218
| 
219
{| style="text-align: left; margin:auto;width: 100%;" 
220
|-
221
| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0</math>
222
|}
223
| style="width: 5px;text-align: right;white-space: nowrap;" |  (15b)
224
|}
225
226
We will neglect hereonwards the third integral in Eq.(15b) by assuming that <math display="inline">r_{m_i}</math> is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(15a) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are
227
228
{| class="formulaSCP" style="width: 100%; text-align: left;" 
229
|-
230
| 
231
{| style="text-align: left; margin:auto;width: 100%;" 
232
|-
233
| style="text-align: center;" | <math>\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)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega -  \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p  d\Gamma +\qquad \\ \displaystyle + \int _{\Omega } {h_j\over 2}{\partial \delta v_i \over \partial x_j} r_{m_i} d\Omega =0\qquad  \end{array}</math>
234
|}
235
| style="width: 5px;text-align: right;white-space: nowrap;" |  (16a)
236
|}
237
238
{| class="formulaSCP" style="width: 100%; text-align: left;" 
239
|-
240
| 
241
{| style="text-align: left; margin:auto;width: 100%;" 
242
|-
243
| style="text-align: center;" | <math>\int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0</math>
244
|}
245
| style="width: 5px;text-align: right;white-space: nowrap;" |  (16b)
246
|}
247
248
In the derivation of the viscous term in Eq.(16a) we have used the following identity (prior to the integration by parts)
249
250
{| class="formulaSCP" style="width: 100%; text-align: left;" 
251
|-
252
| 
253
{| style="text-align: left; margin:auto;width: 100%;" 
254
|-
255
| style="text-align: center;" | <math>{\partial s_{ij} \over \partial x_j}=2\mu {\partial \varepsilon _{ij} \over \partial x_j}=\mu {\partial ^2v_i\over \partial x_j \partial x_j} </math>
256
|}
257
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
258
|}
259
260
Eq.(17) is identically true for the exact incompressible limit <math display="inline">\left({\partial v_i \over \partial x_i}=0\right)</math>.
261
262
===2.2 Convective and pressure gradient projections===
263
264
The computation of the residual terms can be simplified if we introduce now the convective and pressure gradient projections <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively defined as
265
266
{| class="formulaSCP" style="width: 100%; text-align: left;" 
267
|-
268
| 
269
{| style="text-align: left; margin:auto;width: 100%;" 
270
|-
271
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle c_i = r_{m_i} -\rho v_j {\partial v_i \over \partial x_j}\\ \displaystyle \pi _i = r_{m_i} - {\partial p \over \partial x_i} \end{array} </math>
272
|}
273
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
274
|}
275
276
We can express <math display="inline">r_{m_i}</math> in  Eqs.(16a) and (16b) in terms of <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number of  equations by imposing that the residual <math display="inline">r_{m_i}</math> vanishes (in average sense) for both forms given by Eqs.(18). This gives the final system of governing equation as:
277
278
{| class="formulaSCP" style="width: 100%; text-align: left;" 
279
|-
280
| 
281
{| style="text-align: left; margin:auto;width: 100%;" 
282
|-
283
| style="text-align: center;" | <math>\int _\Omega \left[\delta v_i\rho \left({\partial v_i \over \partial t}+v_j {\partial {v_i} \over \partial x_j}\right)+ {\partial \delta v_i \over \partial x_j}\left(\mu {\partial v_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega -  \int _{\Omega } \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_i^p d\Gamma +</math>
284
|-
285
| style="text-align: center;" | <math> + \int _{\Omega } {h_k\over 2}{\partial (\delta v_i) \over \partial x_k} \left(\rho v_j {\partial {v_i} \over \partial x_j} + c_i\right)d\Omega =0\qquad  </math>
286
|}
287
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
288
|}
289
290
{| class="formulaSCP" style="width: 100%; text-align: left;" 
291
|-
292
| 
293
{| style="text-align: left; margin:auto;width: 100%;" 
294
|-
295
| style="text-align: center;" | <math>\int _\Omega q {\partial v_i \over \partial x_i} d\Omega + \int _\Omega \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0 </math>
296
|}
297
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
298
|}
299
300
{| class="formulaSCP" style="width: 100%; text-align: left;" 
301
|-
302
| 
303
{| style="text-align: left; margin:auto;width: 100%;" 
304
|-
305
| style="text-align: center;" | <math>\int _\Omega \delta c_i \rho \left(\rho v_j {\partial {v_i} \over \partial x_j} + c_i\right)d\Omega =0 \qquad \hbox{no sum in }i </math>
306
|}
307
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
308
|}
309
310
{| class="formulaSCP" style="width: 100%; text-align: left;" 
311
|-
312
| 
313
{| style="text-align: left; margin:auto;width: 100%;" 
314
|-
315
| style="text-align: center;" | <math>\int _\Omega \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0\qquad \hbox{no sum in }i </math>
316
|}
317
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
318
|}
319
320
with <math display="inline">i,j,k=1,n_d</math>.  In Eqs.(21) and (22) <math display="inline">\delta c_i</math> and <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\rho </math> and <math display="inline">\tau _i</math> weights are introduced for convenience.
321
322
We note that accounting for the convective and pressure gradient projections enforces the consistency of the formulation as it ensures that the stabilization terms in Eqs.(19&#8211;22) have a residual form which vanishes for the “exact” solution. Neglecting these terms lowers the accuracy of the numerical solution and it makes the formulation more sensitive to the value of the stabilization parameters as shown in <span id='citeF-34'></span>[[#cite-34|[34,36,37]]].
323
324
==3 FINITE ELEMENT DISCRETIZATION==
325
326
We choose <math display="inline">C^\circ </math> continuous linear interpolations of the velocities, the pressure, the convection projections <math display="inline">c_i</math> and the pressure gradient projections <math display="inline">\pi _i</math> over three node triangles (2D) and four node tetrahedra (3D) <span id='citeF-8'></span>[[#cite-8|[8]]]. The linear interpolations are written as
327
328
{| class="formulaSCP" style="width: 100%; text-align: left;" 
329
|-
330
| 
331
{| style="text-align: left; margin:auto;width: 100%;" 
332
|-
333
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle v_i =  N^k \bar v_i^k \quad , \quad p = N^k \bar p^k\\ \displaystyle c_i = N^k \bar c_i^k \quad , \quad \pi _i =  N^k \bar \pi _i^k \end{array} </math>
334
|}
335
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
336
|}
337
338
where the sum goes over the number of nodes of each element <math display="inline">n</math> (<math display="inline">n=3/4</math> for triangles/tetrahedra), <math display="inline">\bar {(\cdot )}^k</math> denotes nodal variables and <math display="inline">N^k</math> are the linear shape functions <span id='citeF-8'></span>[[#cite-8|[8]]].
339
340
Substituting the approximations (23) into Eqs.(19&#8211;22) and choosing the Galerking form with <math display="inline">\delta v_i =q=\delta c_i=\delta \pi _i =N^i</math> leads to following system of discretized equations
341
342
{| class="formulaSCP" style="width: 100%; text-align: left;" 
343
|-
344
| 
345
{| style="text-align: left; margin:auto;width: 100%;" 
346
|-
347
| style="text-align: center;" | <math>\displaystyle {\bf M}\dot{\bar {\bf v}} + {\bf H} \bar {\bf v} - {\bf G}\bar {\bf p}+{\bf C}\bar {\bf c}={\bf f}</math>
348
|}
349
| style="width: 5px;text-align: right;white-space: nowrap;" | (24a)
350
|}
351
352
{| class="formulaSCP" style="width: 100%; text-align: left;" 
353
|-
354
| 
355
{| style="text-align: left; margin:auto;width: 100%;" 
356
|-
357
| style="text-align: center;" | <math>\displaystyle {\bf G}^T \bar {\bf v} + \hat{\bf L}\bar {\bf p}+{\bf Q}\bar {\boldsymbol \pi}={\bf 0}</math>
358
|}
359
| style="width: 5px;text-align: right;white-space: nowrap;" | (24b)
360
|}
361
362
{| class="formulaSCP" style="width: 100%; text-align: left;" 
363
|-
364
| 
365
{| style="text-align: left; margin:auto;width: 100%;" 
366
|-
367
| style="text-align: center;" | <math>\displaystyle \hat {\bf C}\bar{\bf v}+ {\bf M}\bar {\bf c}={\bf 0}</math>
368
|}
369
| style="width: 5px;text-align: right;white-space: nowrap;" | (24c)
370
|}
371
372
{| class="formulaSCP" style="width: 100%; text-align: left;" 
373
|-
374
| 
375
{| style="text-align: left; margin:auto;width: 100%;" 
376
|-
377
| style="text-align: center;" | <math>\displaystyle {\bf Q}^T \bar {\bf p} + \hat {\bf M}\bar {\boldsymbol \pi}={\bf 0}</math>
378
|}
379
| style="width: 5px;text-align: right;white-space: nowrap;" | (24d)
380
|}
381
382
The form of the different matrices is given in the Appendix.
383
384
The solution in time  of the system of Eqs.(24) can be written in general form as
385
386
{| class="formulaSCP" style="width: 100%; text-align: left;" 
387
|-
388
| 
389
{| style="text-align: left; margin:auto;width: 100%;" 
390
|-
391
| style="text-align: center;" | <math>{\bf M} \displaystyle{1\over \Delta t} (\bar{\bf v}^{n+1}-\bar{\bf v}^{n}) + {\bf H}^{n+\theta}\bar{\bf v} ^{n+\theta} - {\bf G}\bar{\bf p}^{n+\theta}+ {\bf C}^{n+\theta}\bar {\bf c}^{n+\theta}={\bf f}^{n+\theta}</math>
392
|}
393
| style="width: 5px;text-align: right;white-space: nowrap;" |  (25a)
394
|}
395
396
{| class="formulaSCP" style="width: 100%; text-align: left;" 
397
|-
398
| 
399
{| style="text-align: left; margin:auto;width: 100%;" 
400
|-
401
| style="text-align: center;" | <math>{\bf G}^T \bar {\bf v}^{n+\theta}+\hat {\bf L}^{n+\theta} \bar {\bf p}^{n+\theta}+{\bf Q}\bar {\boldsymbol  \pi}^{n+\theta}={\bf 0}</math>
402
|}
403
| style="width: 5px;text-align: right;white-space: nowrap;" |  (25b)
404
|}
405
406
{| class="formulaSCP" style="width: 100%; text-align: left;" 
407
|-
408
| 
409
{| style="text-align: left; margin:auto;width: 100%;" 
410
|-
411
| style="text-align: center;" | <math>\hat {\bf C}^{n+\theta} \bar {\bf v}^{n+\theta}+ {\bf M} \bar {\bf c}^{n+\theta}={\bf 0}</math>
412
|}
413
| style="width: 5px;text-align: right;white-space: nowrap;" |  (25c)
414
|}
415
416
{| class="formulaSCP" style="width: 100%; text-align: left;" 
417
|-
418
| 
419
{| style="text-align: left; margin:auto;width: 100%;" 
420
|-
421
| style="text-align: center;" | <math>{\bf G}^T \bar {\bf p}^{n+\theta}+\hat {\bf M}^{n+\theta}\bar {\boldsymbol \pi}^{n+\theta}={\bf 0}</math>
422
|}
423
| style="width: 5px;text-align: right;white-space: nowrap;" |  (25d)
424
|}
425
426
where <math display="inline">{\bf H}^{n+\theta}={\bf H} (\bar{\bf v}^{n+\theta})</math>, <math display="inline">\bar{\bf v}^{n+\theta}</math> are the velocities evaluated at time <math display="inline">n+\theta </math> and the parameter <math display="inline">\theta \in [0,1]</math>. The direct monolitic solution of Eqs.(25) is possible using an adequate iterative scheme. However, we have found more convenient to  use  a fractional step method as described in the next section.
427
428
===3.1 Fractional step method===
429
430
A fractional step scheme is derived by noting that the discretized momentum equation (25a) can be split into the two following equations
431
432
{| class="formulaSCP" style="width: 100%; text-align: left;" 
433
|-
434
| 
435
{| style="text-align: left; margin:auto;width: 100%;" 
436
|-
437
| style="text-align: center;" | <math>{\bf M} \displaystyle{1\over \Delta t} (\tilde{\bf v}^{n+1}-\bar{\bf v}^{n}) + {\bf H}^{n+\theta}\bar{\bf v} ^{n+\theta} - \alpha {\bf G}\bar{\bf p}^{n} + {\bf C}^{n+\theta}\bar {\bf c}^{n+\theta}={\bf f}^{n+\theta}</math>
438
|}
439
| style="width: 5px;text-align: right;white-space: nowrap;" |  (26)
440
|}
441
442
{| class="formulaSCP" style="width: 100%; text-align: left;" 
443
|-
444
| 
445
{| style="text-align: left; margin:auto;width: 100%;" 
446
|-
447
| style="text-align: center;" | <math>{\bf M} \displaystyle{1\over \Delta t} (\bar{\bf v}^{n+1}-\tilde{\bf v}^{n+1})- {\bf G}(\bar{\bf p}^{n+1}-\alpha \bar {\bf p}^{n})={\bf 0}</math>
448
|}
449
| style="width: 5px;text-align: right;white-space: nowrap;" |  (27)
450
|}
451
452
In above equation <math display="inline">\tilde{v}^{n+1}</math> is a predicted value of the velocity at time <math display="inline">n+1</math> and <math display="inline">\alpha </math> is a variable whose values of interest are zero and one. For <math display="inline">\alpha =0</math> (first order scheme) the splitting error is of order <math display="inline"> 0 (\Delta t)</math>, whereas for <math display="inline">\alpha =1</math> (second order scheme) the error is of order <math display="inline">0 (\Delta t^2)</math> <span id='citeF-45'></span>[[#cite-45|[45]]].
453
454
Eqs.(26) and (27) are completed with the following three equations emanating from Eqs.(25b-d)
455
456
{| class="formulaSCP" style="width: 100%; text-align: left;" 
457
|-
458
| 
459
{| style="text-align: left; margin:auto;width: 100%;" 
460
|-
461
| style="text-align: center;" | <math>{\bf G}^T\bar{\bf v}^{n+1}+\hat {\bf L}^n \bar{\bf p}^{n+1}+{\bf Q}\bar {\boldsymbol \pi}^{n}= {\bf 0}</math>
462
|}
463
| style="width: 5px;text-align: right;white-space: nowrap;" |  (28a)
464
|}
465
466
{| class="formulaSCP" style="width: 100%; text-align: left;" 
467
|-
468
| 
469
{| style="text-align: left; margin:auto;width: 100%;" 
470
|-
471
| style="text-align: center;" | <math>\hat {\bf C}^{n+1} \bar{\bf v}^{n+1} + {\bf M} \bar {\bf c}^{n+1}= {\bf 0}</math>
472
|}
473
| style="width: 5px;text-align: right;white-space: nowrap;" |  (28b)
474
|}
475
476
{| class="formulaSCP" style="width: 100%; text-align: left;" 
477
|-
478
| 
479
{| style="text-align: left; margin:auto;width: 100%;" 
480
|-
481
| style="text-align: center;" | <math>{\bf Q}^T \bar{\bf p}^{n+1}+ \hat {\bf M}^{n+1} \bar {\boldsymbol \pi}^{n+1}= {\bf 0}</math>
482
|}
483
| style="width: 5px;text-align: right;white-space: nowrap;" |  (28c)
484
|}
485
486
The value of <math display="inline">\bar{\bf v}^{n+1}</math> obtained from Eq.(27) is substituted into Eq.(28a) to give
487
488
{| class="formulaSCP" style="width: 100%; text-align: left;" 
489
|-
490
| 
491
{| style="text-align: left; margin:auto;width: 100%;" 
492
|-
493
| style="text-align: center;" | <math>{\bf G}^T\tilde{\bf v}^{n+1}+ \Delta t {\bf G}^T {\bf M}^{-1} {\bf G} (\bar{\bf p}^{n+1} - \alpha \bar{\bf p}^n)+ \hat {\bf L}^n {\bf p}^{n+1}+{\bf Q}\bar {\boldsymbol \pi}^n={\bf 0}</math>
494
|}
495
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
496
|}
497
498
The product <math display="inline">{\bf G}^T {\bf M}^{-1} {\bf G}</math> can be approximated by a laplacian matrix, i.e.
499
500
{| class="formulaSCP" style="width: 100%; text-align: left;" 
501
|-
502
| 
503
{| style="text-align: left; margin:auto;width: 100%;" 
504
|-
505
| style="text-align: center;" | <math>{\bf G}^T {\bf M}^{-1} {\bf G}={1\over \rho } {L}\quad \hbox{with }L^{ab}=\int _{ \Omega ^e} {\boldsymbol \nabla }^T N^a {\boldsymbol \nabla } N^b d\Omega  </math>
506
|}
507
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
508
|}
509
510
where <math display="inline">L^{ab}</math> are the element contributions to <math display="inline">{L}</math> (see Appendix).
511
512
The steps of the fractional step scheme chosen here are:
513
514
'''Step 1'''
515
516
The fractional nodal velocities <math display="inline">\tilde {v}^{n+1}</math> can be  explicitely computed from Eq.(26) by
517
518
{| class="formulaSCP" style="width: 100%; text-align: left;" 
519
|-
520
| 
521
{| style="text-align: left; margin:auto;width: 100%;" 
522
|-
523
| style="text-align: center;" | <math>\tilde{\bf v}^{n+1} = \bar{\bf v}^{n} - \Delta t  {\bf M}^{-1}_d [\tilde{\bf H}^{n}\bar{\bf v}^{n} - \alpha {\bf G} \bar{\bf p}^n + {\bf C}^{n}\bar{\bf c}^{n} - {\bf f}^{n}]</math>
524
|}
525
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
526
|}
527
528
where <math display="inline">{\bf M}_d</math> is the diagonal form of <math display="inline">{\bf M}</math> obtaining by lumping the row terms into the corresponding diagonal terms.
529
530
<br/>
531
532
'''Step 2''' Compute <math display="inline">\bar{\bf p}^{n+1}</math> from Eq.(29) as
533
534
{| class="formulaSCP" style="width: 100%; text-align: left;" 
535
|-
536
| 
537
{| style="text-align: left; margin:auto;width: 100%;" 
538
|-
539
| style="text-align: center;" | <math>\bar{\bf p}^{n+1}= -[\hat{\bf L}^n + {\Delta t \over \rho} {\bf L}]^{-1} [{\bf G}^T\tilde{\bf v}^{n+1} - \alpha {\Delta t \over \rho}{\bf L}\bar{\bf p}^n +{\bf Q} \bar{\boldsymbol \pi }^{n}] </math>
540
|}
541
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
542
|}
543
544
'''Step 3''' Compute <math display="inline"> \bar{\bf v}^{n+1}</math> explicitely from Eq.(28a) as
545
546
{| class="formulaSCP" style="width: 100%; text-align: left;" 
547
|-
548
| 
549
{| style="text-align: left; margin:auto;width: 100%;" 
550
|-
551
| style="text-align: center;" | <math>\bar {\bf v}^{n+1}=\tilde{\bf v}^{n+1}+ \Delta t  {\bf M}_d^{-1} {\bf G} (\bar {\bf p}^{n+1}- \alpha \bar {\bf p}^n) </math>
552
|}
553
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
554
|}
555
556
<br/>
557
558
'''Step 4''' Compute <math display="inline"> \bar{\bf c}^{n+1}</math> explicitely from Eq.(28b) as
559
560
{| class="formulaSCP" style="width: 100%; text-align: left;" 
561
|-
562
| 
563
{| style="text-align: left; margin:auto;width: 100%;" 
564
|-
565
| style="text-align: center;" | <math>\bar{\bf c}^{n+1}=- {\bf M}_d^{-1}\hat {\bf C}^{n+1}\bar{\bf v}^{n+1} </math>
566
|}
567
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
568
|}
569
570
'''Step 5''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.(28c) as
571
572
{| class="formulaSCP" style="width: 100%; text-align: left;" 
573
|-
574
| 
575
{| style="text-align: left; margin:auto;width: 100%;" 
576
|-
577
| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{n+1}=- \hat {\bf M}_d^{-1} {\bf Q}^T \bar {\bf p}^{n+1} </math>
578
|}
579
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
580
|}
581
582
where <math display="inline">\hat {\bf M}_d</math> is the lumped form of <math display="inline">\hat {\bf M}</math>. A standard diagonal lumping procedure based in summing up the terms of each row has been used.
583
584
Note that all steps can be solved explicitely except for the computation of the pressure in Eq.(32) which requires to invert the sum of two laplacian matrices. This can be effectively performed using an iterative solution scheme such as the conjugate gradient method.
585
586
Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix <math display="inline">\hat {\bf L}</math> in Eq.(32). This matrix emanates from the FIC formulation.
587
588
The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">\tilde{\bf v}^{n+1}</math> in Eq.(31). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\bf v}^{n+1}</math> in  step 3. The prescribed pressures at the boundary are imposed by making <math display="inline">\bar{\bf p}^n</math> equal to the prescribed pressure values when solving Eq.(32).
589
590
===3.2 Stokes flow===
591
592
Many metal forming processes can be simulated with the assumption that the convective terms are negligible (Stokes flow) <span id='citeF-8'></span>[[#cite-8|[8,43,44]]]. The Stokes  FIC formulation can be readily obtained simply by neglecting the convective terms in the  Navier-Stokes formulation presented earlier. This also implies neglecting the convective stabilization terms in the momentum equations and, consequently, the convective projection variables are not larger necessary. Also the intrinsic time parameters <math display="inline">\tau _i</math> take now the simpler form (see Eq.(12)):
593
594
{| class="formulaSCP" style="width: 100%; text-align: left;" 
595
|-
596
| 
597
{| style="text-align: left; margin:auto;width: 100%;" 
598
|-
599
| style="text-align: center;" | <math>\tau _i={3h_i^2\over 8\mu } </math>
600
|}
601
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
602
|}
603
604
We note again that in metal forming problems the viscosity <math display="inline">\mu </math> will typically be a function of the strain rate and the yield stress of the material <span id='citeF-8'></span>[[#cite-8|8]].
605
606
The resulting discretized system of equations can be written as (see Eqs.(28))
607
608
{| class="formulaSCP" style="width: 100%; text-align: left;" 
609
|-
610
| 
611
{| style="text-align: left; margin:auto;width: 100%;" 
612
|-
613
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {\bf M}\dot{\bar {\bf v}} + {\bf K}\bar{\bf v} - {\bf G}\bar {\bf p}={\bf f}\\  \displaystyle {\bf G}^T \bar {\bf v} + \hat{\bf L}\bar {\bf p}+{\bf Q}\bar {\boldsymbol \pi }={\bf 0}\\ \displaystyle {\bf Q}^T \bar {\bf p} + \hat {\bf M}\bar {\boldsymbol \pi }={\bf 0} \end{array} </math>
614
|}
615
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
616
|}
617
618
The algorithm of previous section can now be implemented.
619
620
The steady-state form of Eqs.(37) can be expressed in matrix form as
621
622
{| class="formulaSCP" style="width: 100%; text-align: left;" 
623
|-
624
| 
625
{| style="text-align: left; margin:auto;width: 100%;" 
626
|-
627
| style="text-align: center;" | <math>\left[\begin{matrix}{\bf K}&-{\bf G}&{\bf 0}\\ -{\bf G}^T & - \hat{\bf L} &-{\bf Q}\\ {\bf 0}& -{\bf Q}^T&-\hat {\bf M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {\bf v}\\ \bar {\bf p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{\bf f}\\ {\bf 0}\\{\bf 0}\\ \end{matrix}\right\} </math>
628
|}
629
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
630
|}
631
632
The system is symmetric and always positive definite and therefore leads to a non singular solution. This property holds for ''any interpolation function'' chosen for <math display="inline">\bar {\bf v},\bar {\bf p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions <span id='citeF-8'></span>[[#cite-8|[8]]].
633
634
A reduced velocity-pressure formulation can be obtained by elliminating the <math display="inline">\bar \pi </math> variables from the last row of Eq.(38)  <span id='citeF-37'></span>[[#cite-37|[37]]].
635
636
The FIC formulation for Stokes flows is applicable for analysis of quasi/fully incompressible solids. An analogous formulation based on solid mechanics concepts is derived in Section 6.
637
638
==4 THERMAL-COUPLED FLOWS. TREATMENT OF THE FREE BOUNDARY==
639
640
===4.1 Thermal coupled flow===
641
642
The effect of temperature can be easily introduced by solving the equation for heat transport coupled to the fluid flow equations.
643
644
The equation of balance of heat  is written in the FIC formulation as <span id='citeF-30'></span>[[#cite-30|[30]]]
645
646
{| class="formulaSCP" style="width: 100%; text-align: left;" 
647
|-
648
| 
649
{| style="text-align: left; margin:auto;width: 100%;" 
650
|-
651
| style="text-align: center;" | <math>r_\phi - {h_j\over 2} {\partial r_\phi  \over \partial x_j}=0\qquad ,\quad j=1,n_d</math>
652
|}
653
| style="width: 5px;text-align: right;white-space: nowrap;" | (39a)
654
|}
655
656
with
657
658
{| class="formulaSCP" style="width: 100%; text-align: left;" 
659
|-
660
| 
661
{| style="text-align: left; margin:auto;width: 100%;" 
662
|-
663
| style="text-align: center;" | <math>r_\phi :=-\rho c \left({\partial \phi  \over \partial t}+v_i {\partial \phi  \over \partial x_i}\right)+ {\partial  \over \partial x_k} \left(k{\partial \phi  \over \partial x_k}\right)+Q</math>
664
|}
665
| style="width: 5px;text-align: right;white-space: nowrap;" | (39b)
666
|}
667
668
In above <math display="inline">\phi </math> is the temperature, <math display="inline">c</math> and <math display="inline">k</math> are the specific heat and the thermal conductivity of the material, respectively, <math display="inline">Q</math> is the heat source and <math display="inline">h_j</math> are the characteristic length distances which are typical of the FIC formulation <span id='citeF-30'></span>[[#cite-30|[30]]].
669
670
Eq.(39a) is completed with the Dirichlet and Neuman boundary conditions for the heat problem. For details see <span id='citeF-35'></span>[[#cite-35|[35]]].
671
672
The convective velocities <math display="inline">v_i</math> in Eq.(39b) are provided by the solution of the fluid flow problem. As usual in metal forming processes, the heat source <math display="inline">Q</math> is a function of the mechanical work generated in the flow of the material during the forming process. The temperature field affects in turn the flow viscosity via its dependence with the yield stress which is very sensitive to the temperature changes. The solution of the heat transfer equation is therefore fully coupled with that of the fluid flow problem <span id='citeF-30'></span>[[#cite-30|[30,46]]].
673
674
===4.2 Treatment of the free boundary motion===
675
676
For the treatment of the free boundary we use a standard volume of fluid (VOF) technique, also known as pseudo-concentration method or level set technique <span id='citeF-22'></span>[[#cite-22|[22,47,48]]]. In the VOF method the motion of the free boudary is followed by solving the following transport equation (written using the FIC formulation)
677
678
{| class="formulaSCP" style="width: 100%; text-align: left;" 
679
|-
680
| 
681
{| style="text-align: left; margin:auto;width: 100%;" 
682
|-
683
| style="text-align: center;" | <math>r_\beta - \underline{{\bar h_j\over 2} {\partial r_\beta  \over \partial x_j}}{\bar h_j\over 2} {\partial r_\beta  \over \partial x_j}=0 \quad \hbox{in }\Omega </math>
684
|}
685
| style="width: 5px;text-align: right;white-space: nowrap;" | (40a)
686
|}
687
688
where
689
690
{| class="formulaSCP" style="width: 100%; text-align: left;" 
691
|-
692
| 
693
{| style="text-align: left; margin:auto;width: 100%;" 
694
|-
695
| style="text-align: center;" | <math>r_\beta :={\partial \beta  \over \partial t}+v_k {\partial \beta  \over \partial x_k} </math>
696
|}
697
| style="width: 5px;text-align: right;white-space: nowrap;" | (40b)
698
|}
699
700
In above <math display="inline">\beta </math> is an auxiliary variable which takes a value equal to one on the free surface and zero elsewhere. The solution of Eq.(40a) in time in the discretized fluid-air domain allows to track the free surface which is characterized by the contours of <math display="inline">\beta </math> which have a unit value.
701
702
The underlined term in Eq.(40a) emerges from the FIC formulation <span id='citeF-33'></span>[[#cite-33|[33]]]. Again, the <math display="inline">\bar h_j</math>'s in Eq.(40a) are characteristic length distances of the order of the element size. This term introduces the necessary stabilization in the transient solution of Eq.(40a).
703
704
===4.3 Fully coupled algorithm for mould filling problems===
705
706
The analysis of mould filling problems typically requires the solution of the coupled-thermal flow problem accounting for the transport of the free surface.
707
708
The simplest explicit algorithm is as follows:
709
710
*  For each time step:
711
* Compute the velocities and the pressure in the fluid <math display="inline">(\bar{v}^{n+1}_i, \bar{p}^{n+1})</math> using the fractional step scheme of Section 3.1.
712
* Compute the nodal temperatures <math display="inline">(\phi ^{n+1})</math> solving Eq.(39a).
713
* Compute the new position of the free surface <math display="inline">(\beta ^{n+1})</math> solving Eq.(40a).
714
715
A number of implicit versions of the algorithm are possible and they all involve an iteration loop within each time step until convergence for the flow variables, the temperature and the free surface position is found. For details see <span id='citeF-47'></span>[[#cite-47|[47,48]]].
716
717
==5 COMPUTATION OF THE CHARACTERISTIC LENGTHS==
718
719
The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Excellent results have been obtained in all problems solved  using linear  tetrahedra with the characteristic length vector  defined by
720
721
{| class="formulaSCP" style="width: 100%; text-align: left;" 
722
|-
723
| 
724
{| style="text-align: left; margin:auto;width: 100%;" 
725
|-
726
| style="text-align: center;" | <math>{\bf h}=h_s {{\bf v}\over {v}}+h_{c} {{\boldsymbol \nabla } v\over \vert{\boldsymbol \nabla }v\vert }</math>
727
|}
728
| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)
729
|}
730
731
where <math display="inline">v=\vert {\bf v}\vert </math> and <math display="inline">h_s</math> and <math display="inline">h_{c}</math> are the “streamline” and “cross wind” contributions given by
732
733
{| class="formulaSCP" style="width: 100%; text-align: left;" 
734
|-
735
| 
736
{| style="text-align: left; margin:auto;width: 100%;" 
737
|-
738
| style="text-align: center;" | <math>h_s=\max ({l}^T_j {\bf v})/{v} </math>
739
|}
740
| style="width: 5px;text-align: right;white-space: nowrap;" | (41b)
741
|}
742
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>h_{c}=\max ({\bf l}^T_j {\boldsymbol \nabla }v)/ \vert {\boldsymbol \nabla }v\vert \quad , \quad  j=1,n_s</math>
749
|}
750
| style="width: 5px;text-align: right;white-space: nowrap;" | (41c)
751
|}
752
753
where <math display="inline">{\bf l}_j</math> are the vectors defining the element sides (<math display="inline">n_s=6</math> for tetrahedra).
754
755
As for the free surface equation the following value of the characteristic length vector <math display="inline">\bar{\bf h}</math> in Eq.(40a) has been taken
756
757
{| class="formulaSCP" style="width: 100%; text-align: left;" 
758
|-
759
| 
760
{| style="text-align: left; margin:auto;width: 100%;" 
761
|-
762
| style="text-align: center;" | <math>\bar{\bf h} =\bar h_s {{\bf v}\over {v}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert } </math>
763
|}
764
| style="width: 5px;text-align: right;white-space: nowrap;" | (42a)
765
|}
766
767
The streamline parameter <math display="inline">\bar h_s</math> has been obtained by Eq.(41b) whereas the cross wind parameter <math display="inline">\bar h_c</math> has been computed by
768
769
{| class="formulaSCP" style="width: 100%; text-align: left;" 
770
|-
771
| 
772
{| style="text-align: left; margin:auto;width: 100%;" 
773
|-
774
| style="text-align: center;" | <math>\bar h_c = \max [{\bf l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3</math>
775
|}
776
| style="width: 5px;text-align: right;white-space: nowrap;" | (42b)
777
|}
778
779
The cross-wind terms in eqs.(41a) and (42a) account for the effect of the gradient of the solution in the stabilization parameters. This is a standard assumption in most “shock-capturing” stabilization procedures <span id='citeF-8'></span>[[#cite-8|[8,24]]].
780
781
==6 FIC FORMULATION FOR QUASI/FULLY INCOMPRESSIBLE  SOLIDS==
782
783
===6.1 Equilibrium equations===
784
785
As usual the governing equations of solid mechanics are written in a Lagrangian reference frame. Following the arguments of   Section 2 the equilibrium equations  for a  solid are written using the FIC technique as <span id='citeF-37'></span>[[#cite-37|[37]]]
786
787
<span id="eq-43"></span>
788
{| class="formulaSCP" style="width: 100%; text-align: left;" 
789
|-
790
| 
791
{| style="text-align: left; margin:auto;width: 100%;" 
792
|-
793
| style="text-align: center;" | <math>{r}_{i}- \frac{{h}_k}{2}{\partial {r}_{i} \over \partial {x}_k}=0 \quad \mbox{in} \;\,\Omega \qquad k=1,n_d </math>
794
|}
795
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
796
|}
797
798
where for the dynamic case
799
800
<span id="eq-44"></span>
801
{| class="formulaSCP" style="width: 100%; text-align: left;" 
802
|-
803
| 
804
{| style="text-align: left; margin:auto;width: 100%;" 
805
|-
806
| style="text-align: center;" | <math>{r}_{i}: =-\rho {\partial ^2 u_i\over \partial t^2} + {\partial {\sigma }_{ij} \over \partial {x}_j}+{b}_i\qquad , \quad j=1,n_d </math>
807
|}
808
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
809
|}
810
811
In  Eqs.([[#eq-43|43]]) and ([[#eq-44|44]]) <math display="inline">u_i</math> are the displacements, <math>{\sigma }_{ij}</math>and <math>{b}_i</math>are the stresses and the body forces, respectively and <math>{h}_k</math>are characteristic length distances  of an arbitrary prismatic domain where equilibrium of forces is considered.
812
813
Equations ([[#eq-43|43]]) and ([[#eq-44|44]]) are completed with the boundary conditions on the displacements <math display="inline">u_i</math>
814
815
<span id="eq-45"></span>
816
{| class="formulaSCP" style="width: 100%; text-align: left;" 
817
|-
818
| 
819
{| style="text-align: left; margin:auto;width: 100%;" 
820
|-
821
| style="text-align: center;" | <math>{u}_i-{u}_i^p =0\quad \mbox{on} \;\,\Gamma _u </math>
822
|}
823
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
824
|}
825
826
and the equilibrium of surface tractions
827
828
<span id="eq-46"></span>
829
{| class="formulaSCP" style="width: 100%; text-align: left;" 
830
|-
831
| 
832
{| style="text-align: left; margin:auto;width: 100%;" 
833
|-
834
| style="text-align: center;" | <math>{\sigma }_{ij}{n}_j- t_i^p - \underline{\frac{1}{2}{h}_k{n}_k{r}_{i}}\frac{1}{2}{h}_k{n}_k{r}_{i} =0\quad \mbox{on} \;\,\Gamma _t </math>
835
|}
836
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
837
|}
838
839
In the above <math display="inline">{u}_i^p</math> and <math display="inline">t_i</math> are prescribed displacements and tractions over the boundaries <math display="inline">\Gamma _u</math>and <math display="inline">\Gamma _t</math>, respectively and <math display="inline">n_i</math> are the components of the unit normal vector.
840
841
Note that for consistency with the  fluid flow equations of previous section and differently from the tradition in solid mechanics, the compression  pressure has been taken as positive.
842
843
===6.2 Stabilized pressure constitutive equations===
844
845
For simplicity the treatment of the constitutive equations  will be explained for the linear elastic model. The approach extends naturally to the non linear elasto-plastic/viscoplastic constitutive equations typical of metal forming problems.
846
847
As usual in quasi-incompressible problems the stresses are split into deviatoric and volumetric (pressure) parts
848
849
<span id="eq-47"></span>
850
{| class="formulaSCP" style="width: 100%; text-align: left;" 
851
|-
852
| 
853
{| style="text-align: left; margin:auto;width: 100%;" 
854
|-
855
| style="text-align: center;" | <math>{\sigma }_{ij}= s_{ij} -p{\delta }_{ij} </math>
856
|}
857
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
858
|}
859
860
where <math>{\delta }_{ij}</math>is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses <math display="inline">s_{ij}</math> are written as
861
862
<span id="eq-48"></span>
863
{| class="formulaSCP" style="width: 100%; text-align: left;" 
864
|-
865
| 
866
{| style="text-align: left; margin:auto;width: 100%;" 
867
|-
868
| style="text-align: center;" | <math>s_{ij}\!\!=\!\! 2G \left({\varepsilon }_{ij}-\frac{1}{3}{\varepsilon }_{v}{\delta }_{ij}\right) </math>
869
|}
870
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
871
|}
872
873
where <math display="inline">G</math>is the shear modulus,
874
875
{| class="formulaSCP" style="width: 100%; text-align: left;" 
876
|-
877
| 
878
{| style="text-align: left; margin:auto;width: 100%;" 
879
|-
880
| style="text-align: center;" | <math>{\varepsilon }_{ij}= \frac{1}{2} \left( {\partial {u}_i \over \partial {x}_j}+ {\partial {u}_j \over \partial x_i} \right) \quad  \mbox{and} \quad  {\varepsilon }_{v}={\varepsilon }_{ii}\,. </math>
881
|}
882
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
883
|}
884
885
The constitutive equation for the pressure <math display="inline">p</math> can be written  for an arbitrary domain of finite size  of volume <math display="inline">V</math> as
886
887
<span id="eq-50"></span>
888
{| class="formulaSCP" style="width: 100%; text-align: left;" 
889
|-
890
| 
891
{| style="text-align: left; margin:auto;width: 100%;" 
892
|-
893
| style="text-align: center;" | <math>{1\over K} p^{av}= - {\Delta V\over V}  </math>
894
|}
895
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
896
|}
897
898
where <math display="inline">K</math> is the bulk modulus of the material and <math display="inline">p^{av}</math> is the average value of the pressure over  domain <math display="inline">V</math>.
899
900
The value of <math display="inline">p^{av}</math> can be approximated as <span id='citeF-37'></span>[[#cite-37|[37]]]
901
902
<span id="eq-51"></span>
903
{| class="formulaSCP" style="width: 100%; text-align: left;" 
904
|-
905
| 
906
{| style="text-align: left; margin:auto;width: 100%;" 
907
|-
908
| style="text-align: center;" | <math>p^{av}: ={1\over V} \int _V p\,dV =p - {h_k\over 2}{\partial p \over \partial x_k} + O (h_k)^2\quad ,\quad  k=1,n_d  </math>
909
|}
910
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
911
|}
912
913
where <math display="inline">p</math> is the pressure at an arbitrary point within the domain <math display="inline">V</math> and <math display="inline">h_k</math> are characteristic lengths of such a domain.
914
915
The ratio <math display="inline">{\Delta V\over V}</math> can be expressed as
916
917
<span id="eq-52"></span>
918
{| class="formulaSCP" style="width: 100%; text-align: left;" 
919
|-
920
| 
921
{| style="text-align: left; margin:auto;width: 100%;" 
922
|-
923
| style="text-align: center;" | <math>{\Delta V\over V} = {\varepsilon }_{v}- {h_k\over 2}{\partial {\varepsilon }_{v} \over \partial x_k} + O (h_k)^2\quad ,\quad k=1,n_d  </math>
924
|}
925
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
926
|}
927
928
Substituting Eqs.([[#eq-51|51]]) and ([[#eq-52|52]]) into Eq.([[#eq-50|50]]) and neglecting second order terms in <math display="inline">h_k</math> gives the FIC constitutive equation for the pressure as
929
930
<span id="eq-53"></span>
931
{| class="formulaSCP" style="width: 100%; text-align: left;" 
932
|-
933
| 
934
{| style="text-align: left; margin:auto;width: 100%;" 
935
|-
936
| style="text-align: center;" | <math>\left({p\over K} + {\varepsilon }_{v}\right)- \frac{{h}_k}{2}{\partial  \over \partial {x}_k} \left({p\over K}  + {\varepsilon }_{v}\right)=0\, ,\quad k=1,n_d  </math>
937
|}
938
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
939
|}
940
941
Note that for <math display="inline">h_k \to 0</math> the standard relationship between the pressure and the volumetric strain of the infinitesimal theory <math display="inline">(p=-K\varepsilon _v )</math> is found.
942
943
For an incompressible material <math display="inline">K\to \infty </math> and Eq.([[#eq-53|53]]) yields
944
945
{| class="formulaSCP" style="width: 100%; text-align: left;" 
946
|-
947
| 
948
{| style="text-align: left; margin:auto;width: 100%;" 
949
|-
950
| style="text-align: center;" | <math>\varepsilon _v - \frac{{h}_k}{2}{\partial {\varepsilon }_{v} \over \partial {x}_k}=0 </math>
951
|}
952
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
953
|}
954
955
Eq.(54) expresses the limit incompressible behaviour of the solid. This equation is identical to Eq.(2) for incompressible flow problems and there arises from the mass continuity condition <span id='citeF-30'></span>[[#cite-30|[30,32]]].
956
957
By combining Eqs. ([[#eq-43|43]]),  ([[#eq-44|44]]),  ([[#eq-47|47]]), ([[#eq-48|48]]) and ([[#eq-53|53]]) a mixed displacement&#8211;pressure formulation can be written as
958
959
<span id="eq-55"></span>
960
<span id="eq-55"></span>
961
{| class="formulaSCP" style="width: 100%; text-align: left;" 
962
|-
963
| 
964
{| style="text-align: left; margin:auto;width: 100%;" 
965
|-
966
| style="text-align: center;" | <math>-\rho {\partial ^2 u_i\over \partial t^2} +  {\partial \sigma _{ij} \over \partial x_i}+{b}_i-\frac{{h}_k}{2} {\partial r_i \over \partial {x}_k}=0  </math>
967
|}
968
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
969
|}
970
971
<span id="eq-56"></span>
972
{| class="formulaSCP" style="width: 100%; text-align: left;" 
973
|-
974
| 
975
{| style="text-align: left; margin:auto;width: 100%;" 
976
|-
977
| style="text-align: center;" | <math>\left(\frac{p}{K}+{\varepsilon }_{v}\right)-\frac{{h}_k}{2}{\partial  \over \partial {x}_k}\left(\frac{p}{K}+{\varepsilon }_{v}\right)=0 </math>
978
|}
979
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
980
|}
981
982
From the observation of Eq.([[#eq-55|55]]) we can obtain after some algebra <span id='citeF-37'></span>[[#cite-37|[37]]]
983
984
{| class="formulaSCP" style="width: 100%; text-align: left;" 
985
|-
986
| 
987
{| style="text-align: left; margin:auto;width: 100%;" 
988
|-
989
| style="text-align: center;" | <math>{\partial  \over \partial x_k}\left(\frac{p}{K}+{\varepsilon }_{v}\right)\simeq {3h_j\over 8G} {\partial r_k \over \partial x_j} </math>
990
|}
991
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
992
|}
993
994
Substituting Eq.(57) into (53) the mass balance equation can be written as
995
996
{| class="formulaSCP" style="width: 100%; text-align: left;" 
997
|-
998
| 
999
{| style="text-align: left; margin:auto;width: 100%;" 
1000
|-
1001
| style="text-align: center;" | <math>\left(\frac{p}{K}+{\varepsilon }_{v}\right)- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0 \quad \hbox{with }\quad  \tau _i = \displaystyle{3h_i^2\over 8G} </math>
1002
|}
1003
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
1004
|}
1005
1006
In the derivation of Eq.(58) we have neglected the terms involving products <math display="inline">h_ih_j</math> for <math display="inline">h_i \not = h_j</math>.
1007
1008
The pressure gradient projections <math display="inline">\pi _i</math> are introduced now as
1009
1010
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1011
|-
1012
| 
1013
{| style="text-align: left; margin:auto;width: 100%;" 
1014
|-
1015
| style="text-align: center;" | <math>\pi _i \equiv r_i - {\partial p \over \partial x_i} </math>
1016
|}
1017
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
1018
|}
1019
1020
The final system of governing equations is
1021
1022
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1023
|-
1024
| 
1025
{| style="text-align: left; margin:auto;width: 100%;" 
1026
|-
1027
| style="text-align: center;" | <math>\rho {\partial ^2 u_i\over \partial t^2} - {\partial \sigma _{ij} \over \partial x_i} -{b}_i=0 </math>
1028
|}
1029
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
1030
|}
1031
1032
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1033
|-
1034
| 
1035
{| style="text-align: left; margin:auto;width: 100%;" 
1036
|-
1037
| style="text-align: center;" | <math>{\Delta p\over k}+ {\partial (\Delta u_i) \over \partial x_i} + \sum \limits _{i=1}^{n_d} \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)=0 </math>
1038
|}
1039
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
1040
|}
1041
1042
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1043
|-
1044
| 
1045
{| style="text-align: left; margin:auto;width: 100%;" 
1046
|-
1047
| style="text-align: center;" | <math>r_i \equiv {\partial p \over \partial x_i}+\pi _i =0 </math>
1048
|}
1049
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
1050
|}
1051
1052
In Eqs.(60)&#8211;(62) we note the following:
1053
1054
* The stabilization terms have been neglected in Eq.(60). These terms were useful to derive the stabilized form of Eq.(61). However they are not longer necessary as the convective terms do not appear in  the equilibrium equations.
1055
* The constitutive equation for the pressure (Eq.(61)) has been written  in an incremental form. This is more convenient for non linear material behaviour typical of metal forming situations. Here, appropriate elasto-plastic or elasto-viscoplastic constitutive laws relating stresses and strains in incremental or rate forms must be used <span id='citeF-1'></span>[[#cite-1|[1-5,43,44]]].
1056
1057
===6.3 Finite element equations===
1058
1059
The weighted residual form of the governing equations can be written as (after integration by parts of the relevant terms)
1060
1061
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1062
|-
1063
| 
1064
{| style="text-align: left; margin:auto;width: 100%;" 
1065
|-
1066
| style="text-align: center;" | <math>\int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}\,d\Omega + \int _{\Omega }\delta {\varepsilon }_{ij}{\sigma }_{ij}\,d\Omega{-} \int _{\Omega }\delta u_i b_i \,d\Omega  - \int _{\Gamma _t}\delta {u}_it_i^p \,d\Gamma _t=0</math>
1067
|}
1068
| style="width: 5px;text-align: right;white-space: nowrap;" |  (63a)
1069
|}
1070
1071
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1072
|-
1073
| 
1074
{| style="text-align: left; margin:auto;width: 100%;" 
1075
|-
1076
| style="text-align: center;" | <math>\int _{\Omega }q \left({\Delta p\over K} + {\partial (\Delta u_i) \over \partial x_i} \right)\,d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}{\partial q \over \partial x_i} \tau _{i}\left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0</math>
1077
|}
1078
| style="width: 5px;text-align: right;white-space: nowrap;" |  (63b)
1079
|}
1080
1081
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1082
|-
1083
| 
1084
{| style="text-align: left; margin:auto;width: 100%;" 
1085
|-
1086
| style="text-align: center;" | <math>\int _{\Omega }\delta \pi _i \tau _{i} \left({\partial p \over \partial x_i}+\pi _i \right)d\Omega =0\quad \hbox{no sum in }i </math>
1087
|}
1088
| style="width: 5px;text-align: right;white-space: nowrap;" |  (63c)
1089
|}
1090
1091
where again the <math display="inline">\tau _i</math> are introduced in Eq.(63c) for symmetry reasons.
1092
1093
The finite element discretization of the displacements, the pressure and the pressure gradient projections is written  by expressions identical  to Eqs.(23) with the nodal variables now being a function of the time <math display="inline">t</math>. Substituting the approximations into eqs.(63) and using the Galerkin form gives the following system of discretized equations
1094
1095
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1096
|-
1097
| 
1098
{| style="text-align: left; margin:auto;width: 100%;" 
1099
|-
1100
| style="text-align: center;" | <math>{\bf M} \ddot{\bar {\bf u}} + {\bf g} -  {\bf f}={\bf 0}</math>
1101
|}
1102
| style="width: 5px;text-align: right;white-space: nowrap;" | (64a)
1103
|}
1104
1105
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1106
|-
1107
| 
1108
{| style="text-align: left; margin:auto;width: 100%;" 
1109
|-
1110
| style="text-align: center;" | <math>{\bf G}^T \Delta \bar {\bf u} +{\bf C} \Delta \bar {\bf p}+ {\bf L} \bar {\bf p} + {\bf Q}\bar {\boldsymbol \pi }=  {\bf 0}</math>
1111
|}
1112
| style="width: 5px;text-align: right;white-space: nowrap;" |  (64b)
1113
|}
1114
1115
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1116
|-
1117
| 
1118
{| style="text-align: left; margin:auto;width: 100%;" 
1119
|-
1120
| style="text-align: center;" | <math>{\bf Q}^T \bar {\bf p}+ \bar {\bf C} \bar {\boldsymbol \pi }=  {\bf 0}</math>
1121
|}
1122
| style="width: 5px;text-align: right;white-space: nowrap;" |  (64c)
1123
|}
1124
1125
where <math display="inline">\ddot{\bar {u}}</math> is the nodal acceleration vector,
1126
1127
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1128
|-
1129
| 
1130
{| style="text-align: left; margin:auto;width: 100%;" 
1131
|-
1132
| style="text-align: center;" | <math>M_{ij}=\int _{\Omega ^e} \rho N_i N_j\,d\Omega  </math>
1133
|}
1134
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
1135
|}
1136
1137
is the mass matrix,
1138
1139
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1140
|-
1141
| 
1142
{| style="text-align: left; margin:auto;width: 100%;" 
1143
|-
1144
| style="text-align: center;" | <math>{\textbf g}= \int _\Omega {B}^T {\boldsymbol \sigma } d\Omega  </math>
1145
|}
1146
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
1147
|}
1148
1149
is the internal nodal force vector and the rest of matrices and vectors are defined in the Appendix. Note that the expression of <math display="inline">{g}</math> of eq.(66) is adequate for non linear  analysis.
1150
1151
A four steps semi-explicit time integration algorithm can be derived as follows:
1152
1153
1154
'''Step 1.''' Compute the nodal velocities <math display="inline">\dot{\bar {\textbf u}}^{n+1/2}</math> from Eq.(64a)
1155
1156
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1157
|-
1158
| 
1159
{| style="text-align: left; margin:auto;width: 100%;" 
1160
|-
1161
| style="text-align: center;" | <math> \dot{\bar {\textbf u}}^{n+1/2} = \dot{\bar {\textbf u}}^{n-1/2}+\Delta t {\textbf M}^{-1}_d ({\textbf f}^n - {\textbf g}^n)</math>
1162
|}
1163
| style="width: 5px;text-align: right;white-space: nowrap;" |  (67a)
1164
|}
1165
1166
'''Step 2'''. Compute the nodal displacements <math display="inline">\bar {\textbf u}^{n+1}</math>
1167
1168
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1169
|-
1170
| 
1171
{| style="text-align: left; margin:auto;width: 100%;" 
1172
|-
1173
| style="text-align: center;" | <math>\bar {\textbf u}^{n+1} = \bar{\textbf u}^n + \Delta t \dot{\bar {\textbf u}}^{n+1/2} </math>
1174
|}
1175
| style="width: 5px;text-align: right;white-space: nowrap;" |  (67b)
1176
|}
1177
1178
'''Step 3'''. Compute the nodal pressures <math display="inline">\bar {\textbf p}^{n+1}</math> from Eq.(64b)
1179
1180
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1181
|-
1182
| 
1183
{| style="text-align: left; margin:auto;width: 100%;" 
1184
|-
1185
| style="text-align: center;" | <math>\bar {\textbf p}^{n+1}=-[ {\textbf C} + {\textbf L}]^{-1} [\Delta t {\textbf G}^T \dot{\bar {\textbf u}}^{n+1/2}- {\textbf C} \bar {\textbf p}^n + {\textbf Q} \bar {\boldsymbol \pi }^n]</math>
1186
|}
1187
| style="width: 5px;text-align: right;white-space: nowrap;" |  (67c)
1188
|}
1189
1190
'''Step 4'''. Compute the nodal projected pressure gradients <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math> from Eq.(64c)
1191
1192
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1193
|-
1194
| 
1195
{| style="text-align: left; margin:auto;width: 100%;" 
1196
|-
1197
| style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1}=- \bar {\textbf C}_d^{-1} {\textbf Q}^T \bar {\textbf p}^{n+1}</math>
1198
|}
1199
| style="width: 5px;text-align: right;white-space: nowrap;" |  (67d)
1200
|}
1201
1202
In above, all matrices are evaluated at <math display="inline">t^{n+1}</math>, <math display="inline">(\cdot )_d </math>  denotes a lumped diagonal matrix and
1203
1204
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1205
|-
1206
| 
1207
{| style="text-align: left; margin:auto;width: 100%;" 
1208
|-
1209
| style="text-align: center;" | <math>{\textbf g}^n =\int _{\Omega ^e} [{B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
1210
|}
1211
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
1212
|}
1213
1214
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear) constitutive law.
1215
1216
Note that steps 1, 2 and 4 are ''fully explicit'' as a diagonal form of matrices <math display="inline">{\textbf M}</math> and <math display="inline">\bar{\textbf C}</math> has been chosen. The solution of step 3 requires invariably the inverse of a Laplacian matrix. This can be an inexpensive process using an iterative equation solution method (e.g. a preconditioned conjugate gradient method).
1217
1218
For the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">{\textbf C}=0 </math> in all above equations.
1219
1220
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme (see Section 6.5 and <span id='citeF-9'></span>[[#cite-9|[9,37]]]).
1221
1222
'''''Fully explicit algorithm'''''
1223
1224
A fully explicit four steps algorithm can be obtained by computing <math display="inline">\bar {p}^{n+1}</math> from step 3 in eq.(67c) as follows
1225
1226
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1227
|-
1228
| 
1229
{| style="text-align: left; margin:auto;width: 100%;" 
1230
|-
1231
| style="text-align: center;" | <math>\bar{\textbf p}^{n+1}= - {\textbf C}_d^{-1}[\Delta t {\textbf G}^T \dot{\bar {\textbf u}}^{n+1/2} - ({\textbf C}_d - {\textbf L}) \bar{\textbf p}^n + {\textbf Q} \bar{\boldsymbol \pi }^n] </math>
1232
|}
1233
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
1234
|}
1235
1236
Note that the explicit algorithm is not applicable in the full incompressible limit as the solution of Eq.(69) breaks down for <math display="inline">K=\infty </math> and <math display="inline">{\textbf C}=0</math>. The explicit form can however be used with success in problems where quasi-incompressible regions exist adjacent to standard “compressible” zones. Examples of this kind are shown in the next section. In both cases  the semi-explicit and fully explicit schemes gave identical results with important savings in both computer time and memory storage requirements obtained when using the explicit form.
1237
1238
===6.4 Thermal coupled effects===
1239
1240
The effect of temperature can be easily accounted for by solving the heat transfer equation formulated in a Lagrangian frame as
1241
1242
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1243
|-
1244
| 
1245
{| style="text-align: left; margin:auto;width: 100%;" 
1246
|-
1247
| style="text-align: center;" | <math>\rho c {\partial \phi  \over \partial t}-{\partial  \over \partial x_j}\left(k{\partial \phi  \over \partial x_j}\right)- Q=0 </math>
1248
|}
1249
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
1250
|}
1251
1252
As usual, the source <math display="inline">Q</math> is dependent of the mechanical work generated during the forming process.
1253
1254
Note that the convective terms do not enter into Eq.(70) This also eliminates the need to stabilize the numerical solution.
1255
1256
The coupled thermal-mechanical problem requires the computation of the temperature <math display="inline">\phi </math> at each time step using a transient solution scheme <span id='citeF-49'></span>[[#cite-49|[47,48,50]]].
1257
1258
'''Remark'''. Above formulation is similar to that developed by Chiumenti ''et al.'' <span id='citeF-27'></span>[[#cite-27|[27,29]]] and Cervera ''et al.'' <span id='citeF-28'></span>[[#cite-28|[28]]] for analysis of incompressible problems in solid mechanics using a sub-grid scale approach.
1259
1260
===6.5 Computation of the intrinsic time parameter for quasi/fully incompressible solids===
1261
1262
In solid mechanics applications it is usual to accept that all <math display="inline">\tau _i</math> are identical and constant within each element and given by
1263
1264
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1265
|-
1266
| 
1267
{| style="text-align: left; margin:auto;width: 100%;" 
1268
|-
1269
| style="text-align: center;" | <math>\tau ^{(e)} =  {3(h^{(e)})^2\over 8G} \quad \hbox{with}\quad  h^{(e)} =[V^{(e)}]^{1/n_d}\, T </math>
1270
|}
1271
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
1272
|}
1273
1274
where <math display="inline">V^{(e)}</math> is the element volume (or the element area for 2D problems). This expression for <math display="inline">\tau ^{(e)}</math> does not take into account the element distorsions along a particular direction during the deformation process.
1275
1276
The correct value of the shear modulus in the expression of <math display="inline">\tau ^{(e)}</math> is another sensitive issue as, obviously, for non linear problems the value of <math display="inline">G</math> will differ from the elastic modulus. This fact has been identified by Cervera ''et al.'' <span id='citeF-28'></span>[[#cite-28|[28]]] for non linear analysis of incompressible problems using linear triangles.
1277
1278
A useful alternative to compute <math display="inline">\tau ^{(e)}</math> for explicit non linear transient situations is to make use of the value of the speed of sound in an elastic solid, defined by
1279
1280
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1281
|-
1282
| 
1283
{| style="text-align: left; margin:auto;width: 100%;" 
1284
|-
1285
| style="text-align: center;" | <math>c=\sqrt{{E\over \rho }} </math>
1286
|}
1287
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
1288
|}
1289
1290
where <math display="inline">E</math> is the Young modulus. The stability condition for explicit dynamic computations is given by the Courant condition defined as <span id='citeF-8'></span>[[#cite-8|[8]]]
1291
1292
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1293
|-
1294
| 
1295
{| style="text-align: left; margin:auto;width: 100%;" 
1296
|-
1297
| style="text-align: center;" | <math>\Delta t^{(e)}\le \Delta t^{(e)}_c ={h^{(e)}\over c} </math>
1298
|}
1299
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
1300
|}
1301
1302
where <math display="inline">\Delta t^{(e)}_c</math> is the critical time step for the element.
1303
1304
Accepting that <math display="inline">G\simeq {E\over 3}</math> for the incompressible case and using Eqs.(72) and (73) (assuming the identity in Eq.(73)) an alternative expression for the element intrinsic time parameter in terms of the critical time step can be found as
1305
1306
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1307
|-
1308
| 
1309
{| style="text-align: left; margin:auto;width: 100%;" 
1310
|-
1311
| style="text-align: center;" | <math>\tau ^{(e)}={[\Delta t^{(e)}_c]^2\over \rho } </math>
1312
|}
1313
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
1314
|}
1315
1316
Eq.(74) shows clearly that the intrinsic time parameter varies across the mesh as a function of the critical time step for each element.
1317
1318
==7 APPLICATIONS TO CASTING AND BULK METAL FORMING PROCESSES==
1319
1320
===7.1 Aluminium casting simulation===
1321
1322
A numerical simulation of an aluminium casting process is presented as a demonstration of  the  accuracy  of the stabilized formulation. The computations are performed with the finite element code VULCAN where the stabilized FEM presented has been implemented <span id='citeF-56'></span>[[#cite-56|[56]]].
1323
1324
The analysis simulates the casting process of an aluminium (AlSi7Mg) specimen in a steel (X40CrMoV5) mould. Material behavior of aluminium casting has been modeled by a fully coupled thermo-viscoplastic model, while the steel mould  has been modeled by a simpler thermo-elastic model <span id='citeF-51'></span>[[#cite-51|[51]]]. Geometrical and material data were provided by the foundry RUFFINI. Figure 1  shows the finite element mesh used for the part and the cooling system. The full mesh, including the mould has  380.000 four node tetrahedra. The pouring  temperature is 650<math display="inline">^\circ </math>C. Initial temperature for the mould is obtained through a thermal die-cycling simulation. Figure 2 shows the evolution of the mould temperature after 6 cycles. The cooling system has been kept at 20<math display="inline">^\circ </math>C. Filling evolution has been simulated as  in a pressure die-casting process using the stabilized VOF technique described in <span id='citeF-49'></span>[[#cite-49|[47,48]]]. Figure 3 shows different time steps of the simulation.
1325
1326
<div id='img-1'></div>
1327
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1328
|-
1329
|[[File:Draft_Samper_805805629_6310_Fig1cast.png|351px|Finite element discretization of the aluminium casting.]]
1330
|- style="text-align: center; font-size: 75%;"
1331
| colspan="1" | '''Figure 1:''' Finite element discretization of the aluminium casting.
1332
|}
1333
1334
<div id='img-2'></div>
1335
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1336
|-
1337
|[[File:Draft_Samper_805805629_9923_Fig2-cast.png|351px|Temperature die-cycling.]]
1338
|- style="text-align: center; font-size: 75%;"
1339
| colspan="1" | '''Figure 2:''' Temperature die-cycling.
1340
|}
1341
1342
The final temperature field obtained after the filling simulation is taken as the initial condition for the solidification and cooling analysis. Temperature and liquid-fraction distributions during solidification are shown in Figures 4 and 5, respectively. The heat transfer coefficient takes into account the air-gap resistance due to the casting shrinkage during the solidification process. Figure 6 shows volumetric and von Mises deviatoric stress distributions in a x-y section. The figures also show the air-gap between the part and the mould, responsible of a non-uniform heat flux at the contact interface.
1343
1344
Other examples of application of the stabilized formulation to casting problems can be found in <span id='citeF-47'></span>[[#cite-47|[47-51]]].
1345
1346
<div id='img-3'></div>
1347
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1348
|-
1349
|[[Image:Draft_Samper_805805629-MATER.png|300px|]]
1350
|[[Image:Draft_Samper_805805629-MATER-01.png|300px|]]
1351
|-
1352
|[[Image:Draft_Samper_805805629-MATER-02.png|300px|]]
1353
|[[Image:Draft_Samper_805805629-MATER-03.png|300px|]]
1354
|-
1355
| colspan="2"|[[Image:Draft_Samper_805805629-MATER-04.png|300px|Filling evolution: pressure die-casting simulation.]]
1356
|- style="text-align: center; font-size: 75%;"
1357
| colspan="2" | '''Figure 3:''' Filling evolution: pressure die-casting simulation.
1358
|}
1359
1360
<div id='img-4'></div>
1361
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1362
|-
1363
|[[Image:Draft_Samper_805805629-tempe_1.png|210px|]]
1364
|[[Image:Draft_Samper_805805629-tempe_2.png|210px|]]
1365
|[[Image:Draft_Samper_805805629-tempe_3.png|210px|Temperature evolution during cooling phase.]]
1366
|- style="text-align: center; font-size: 75%;"
1367
| colspan="3" | '''Figure 4:''' Temperature evolution during cooling phase.
1368
|}
1369
1370
<div id='img-5'></div>
1371
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1372
|-
1373
|[[Image:Draft_Samper_805805629-solid_frac.png|360px|]]
1374
|[[Image:Draft_Samper_805805629-solid_frac-01.png|360px|]]
1375
|-
1376
|[[Image:Draft_Samper_805805629-solid_frac-02.png|360px|]]
1377
|[[Image:Draft_Samper_805805629-solid_frac-03.png|360px|Liquid-fraction evolution during phase-change.]]
1378
|- style="text-align: center; font-size: 75%;"
1379
| colspan="2" | '''Figure 5:''' Liquid-fraction evolution during phase-change.
1380
|}
1381
1382
<div id='img-6'></div>
1383
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1384
|-
1385
|[[File:Draft_Samper_805805629_9790_Fig6.png|351px|Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).]]
1386
|- style="text-align: center; font-size: 75%;"
1387
| colspan="1" | '''Figure 6:''' Aluminium casting. Stress-trace and von Mises deviatoric stress indicator during  phase-change (plane xy).
1388
|}
1389
1390
===7.2 Side pressing of a cylinder===
1391
1392
A cylinder 100 mm long with a radius of 100 mm is subjected to sidepressing between two plane dies. It is compressed to 100 mm. The material properties are the following: <math display="inline">E =217</math> GPa, <math display="inline">\nu =0.3</math>, <math display="inline">\rho =7830</math> kg/m<math display="inline">^3</math>, <math display="inline">\sigma _0=170</math> MPa, <math display="inline">H=30</math> MPa,  friction coefficient = 0.2. The die velocity is assumed to be 2 m/s. Initial set-up is shown in Figure [[#img-7|7]]. A quarter of a cylinder was discretized with tri-linear hexahedra and linear tetrahedra meshes.
1393
1394
Figure [[#img-8|8]] shows the results obtained using the hexahedral mesh and a standard mixed formulation. The results show the distribution of the effective plastic strain and pressure on the deformed shape. The sensitivity of the FIC results with the expression of the intrinsic time parameter of Eq.(71) was studied by defining <math display="inline">\tau ^{(e)} =a  {[h^{(e)}]^2\over G}</math> and solving the problem for different values of <math display="inline">a</math>. The results obtained with the FIC method are shown in Figs. [[#img-9|9]] and [[#img-10|10]] for <math display="inline">a =0.1</math> and <math display="inline">0.03</math>, respectively. The alternative expression for <math display="inline">\tau ^{(e)}</math> of Eq.(74) has also been studied. The results for this case are shown in Fig. [[#img-11|11]]. Quite a good agreement can be seen between the FIC solutions and the reference solution with the best results for the effective plastic strain obtained for <math display="inline">\tau ^{(e)}</math> calculated according to Eq. (71) with <math display="inline">a=0.03</math>. The results for the two alternative formulae for <math display="inline">\tau ^{(e)}</math> are similar, but those obtained using Eq. (71) seem to be slightly better. In any case, the results on the pressure distribution are quite insentitive to the value of <math display="inline">\tau ^{(e)}</math>. This is also confirmed in Fig. [[#img-12|12]]b, which displays the distribution of the pressure along the line <math display="inline">ABCDEA</math> defined in Fig. [[#img-12|12]]a. A small perturbation can be seen at the sharp edges of the deformed body.
1395
1396
All the calculations here have been carried out using fully explicit version of the algorithm which is more efficient than the semi-implicit one with giving practically the same results.
1397
1398
<div id='img-7'></div>
1399
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1400
|-
1401
|[[Image:Draft_Samper_805805629-cyltn3-pi-v5e5-msh.png|270px|]]
1402
|[[Image:Draft_Samper_805805629-cylhc-msh-new.png|270px|Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh]]
1403
|- style="text-align: center; font-size: 75%;"
1404
| colspan="2" | '''Figure 7:''' Sidepressing of a cylinder: (a) initial tetrahedral mesh; (b) initial hexahedral mesh
1405
|}
1406
1407
<div id='img-8'></div>
1408
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1409
|-
1410
|[[Image:Draft_Samper_805805629-cylhc-ep-new.png|228px|]]
1411
|[[Image:Draft_Samper_805805629-cylhc-hp-new.png|228px|Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution ]]
1412
|- style="text-align: center; font-size: 75%;"
1413
| colspan="2" | '''Figure 8:''' Sidepressing of a cylinder, mixed formulation,  hexahedral mesh: (a) effective plastic strain; (b) pressure distribution 
1414
|}
1415
1416
<div id='img-9'></div>
1417
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1418
|-
1419
|[[Image:Draft_Samper_805805629-cylt-fine-fast-fic01-eps.png|228px|]]
1420
|[[Image:Draft_Samper_805805629-cylt-fine-fast-fic01-p.png|228px|Sidepressing of a cylinder, FIC algorithm  (α=0.1),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
1421
|- style="text-align: center; font-size: 75%;"
1422
| colspan="2" | '''Figure 9:''' Sidepressing of a cylinder, FIC algorithm  (<math>\alpha=0.1</math>),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution 
1423
|}
1424
1425
<div id='img-10'></div>
1426
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1427
|-
1428
|[[Image:Draft_Samper_805805629-cylt-fine-fast-fic-eps.png|228px|]]
1429
|[[Image:Draft_Samper_805805629-cylt-fine-fast-fic-p.png|228px|Sidepressing of a cylinder, FIC algorithm  (a=0.03),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
1430
|- style="text-align: center; font-size: 75%;"
1431
| colspan="2" | '''Figure 10:''' Sidepressing of a cylinder, FIC algorithm  (<math>a=0.03</math>),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution 
1432
|}
1433
1434
<div id='img-11'></div>
1435
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1436
|-
1437
|[[Image:Draft_Samper_805805629-cylt-fine07-fast-v1e12-ab2.png|228px|]]
1438
|[[Image:Draft_Samper_805805629-cylt-fine07-fast-v1e12-p-ab2.png|228px|Sidepressing of a cylinder, FIC algorithm  (τ<sup>(e)</sup> calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution ]]
1439
|- style="text-align: center; font-size: 75%;"
1440
| colspan="2" | '''Figure 11:''' Sidepressing of a cylinder, FIC algorithm  (<math>\tau ^{(e)}</math> calculated according to Eq. (74)),  tetrahedra, mesh  of 22186 elements: (a) effective plastic strain; (b) pressure distribution 
1441
|}
1442
1443
<div id='img-12'></div>
1444
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1445
|-
1446
|[[Image:Draft_Samper_805805629-cylt-linia.png|300px|]]
1447
|[[File:Draft_Samper_805805629_8020_Fig12b.png|350px]]
1448
|- style="text-align: center; font-size: 75%;"
1449
| colspan="2" | '''Figure 12:''' a) Definition of the  line for comparison of pressure distribution b) Pressure distribution along the line <math>ABCDEA</math>
1450
|}
1451
1452
===7.3 Backward extrusion===
1453
1454
Backward extrusion of a cylinder made of steel 16MNCr5 has been analysed using an axi-symmetric formulation. This is a benchmark example of the finite element program for forming simulation MARC/Autoforge <span id='citeF-54'></span>[[#cite-54|[54]]]. The tooling and billet geometry are given in Fig. [[#img-13|13]]a. Initial material dimensions are the following: length 30 mm and diameter 30 mm. The punch of diameter 20 mm has a prescribed stroke of 28 mm. Material properties are as follows: Young's modulus <math display="inline">E=3.24\cdot 10^5</math> MPa, Poissson's coefficient <math display="inline">\nu=0.3</math>, material density <math display="inline">\rho=8120</math> kg/m<math display="inline">^3</math>, yield stress <math display="inline">\sigma _{Y0}=300</math> MPa and hardening modulus <math display="inline">H=50</math> MPa. Friction between the material and tools is defined by the Coulomb friction coefficient <math display="inline">\mu=0.1</math>.
1455
1456
<div id='img-13'></div>
1457
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1458
|-
1459
|[[Image:Draft_Samper_805805629-bext-geo.png|330px|]]
1460
|[[Image:Draft_Samper_805805629-bext4a-quad-nurbs-sprn9-eps.png|266px|]]
1461
|[[Image:Draft_Samper_805805629-bext4a_036-eps.png|279px|]]
1462
|-
1463
|(a)
1464
|colspan="2" |(b)
1465
|- style="text-align: center; font-size: 75%;"
1466
| colspan="3" | '''Figure 13:''' Backward extrusion a) geometry definition. Final deformed shape with effective plastic strain distribution; b)  solution with quadrilaterals and mixed formulation; c)  solution with triangles and the FIC algorithm 
1467
|}
1468
1469
The simulation of the backward extrusion process was carried out with a particularization of the FIC formulation of Section 6 for axisymmetric solids. Regeneration of the  mesh was performed when element distorsion was excessive. Figures [[#img-13|13]]b and [[#img-13|13]]c show the results in the form of the final deformed shape with the distribution of the effective plastic strain obtained using quadrilaterals with a mixed formulation, and using triangles and the FIC algorithm, respectively. The results are in a good agreement with the solution given in <span id='citeF-54'></span>[[#cite-54|[54]]]. This example demonstrates the efficiency of the FIC algorithm for simulation of bulk forming processes. Different stages of the forming process are shown in Fig. [[#img-14|14]].
1470
1471
<div id='img-14'></div>
1472
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1473
|-
1474
|[[Image:Draft_Samper_805805629-bext4a_002-eps.png|330px|]]
1475
|[[Image:Draft_Samper_805805629-bext4a_014-eps.png|330px|]]
1476
|[[Image:Draft_Samper_805805629-bext4a_022-eps.png|330px|]]
1477
|- style="text-align: center; font-size: 75%;"
1478
| colspan="3" | '''Figure 14:''' Backward extrusion &#8211; deformed shapes with effective plastic strain distribution at different stages of forming. Solution with triangles and the FIC algorithm 
1479
|}
1480
1481
===7.4 Forming of a hose-clamp band===
1482
1483
The stabilized formulation was applied to model the manufacturing of a hose-clamp band   of steel AISI 409L (Fig. [[#img-15|15]]). The initial set-up of the tooling and band is shown in Fig. [[#img-16|16]]. A series of grooves are forged in the band by the roll passing over the band placed on the toothed punch. The band thickness is 0.7 mm and its width 8 mm. Plane strain conditions have been assumed. Material properties are as follows: Young's modulus <math display="inline">E=2.1\cdot 10^5</math> MPa, Poissson's coefficient <math display="inline">\nu=0.33</math>, material density <math display="inline">\rho=7800</math> kg/m<math display="inline">^3</math>, the true stress&#8211;true strain relationship is given by the power Ludwik&#8211;Nadai equation <math display="inline">\sigma _{Y}=623(0.36822\cdot 10^{-2}+\varepsilon _p)^{0.1362}</math> MPa, friction between the material and tools is defined by the Coulomb friction coefficient <math display="inline">\mu=0.1</math>.
1484
1485
<div id='img-15'></div>
1486
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1487
|-
1488
|[[Image:Draft_Samper_805805629-hose-clip.png|530px|A hose clamp ]]
1489
|- style="text-align: center; font-size: 75%;"
1490
| colspan="1" | '''Figure 15:''' A hose clamp 
1491
|}
1492
1493
<div id='img-16'></div>
1494
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1495
|-
1496
|[[Image:Draft_Samper_805805629-p1222e_ini.png|410px|A hose clamp &#8211; initial set-up ]]
1497
|- style="text-align: center; font-size: 75%;"
1498
| colspan="1" | '''Figure 16:''' A hose clamp &#8211; initial set-up 
1499
|}
1500
1501
<div id='img-17'></div>
1502
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1503
|-
1504
|[[Image:Draft_Samper_805805629-p1222d_12-eps.png|470px|]]
1505
|-
1506
|[[Image:Draft_Samper_805805629-p1222e_29-eps.png|440px|A hose clamp &#8211; deformed shapes at different stages of forming with distribution of effective plastic strain ]]
1507
|- style="text-align: center; font-size: 75%;"
1508
| '''Figure 17:''' A hose clamp &#8211; deformed shapes at different stages of forming with distribution of effective plastic strain 
1509
|}
1510
1511
Simulation was carried out using linear triangular elements and  the FIC formulation. As in the previous example the  meshes where  regenerated when element distorsion was excessive. The purpose of the simulation was to check if the expected groove depth and tooth height in the band were obtained. Figures [[#img-17|17]]a and [[#img-17|17]]b show the results in the form of the  deformed shape with the distribution of the effective plastic strain at different stages of forming. The results are in good agreement with the real process. Figure  [[#img-18|18]] shows a detail  of the deformed shape with finite element discretization and  distribution of effective plastic strain. In this figure the obtained dimensions are compared with the required ones shown in brackets. Effects of elastic springback are also clearly seen.
1512
1513
<div id='img-18'></div>
1514
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1515
|-
1516
|[[Image:Draft_Samper_805805629-clamp-detail.png|460px|Detail of the deformed shape with finite element discretization and  distribution of effective plastic strain ]]
1517
|- style="text-align: center; font-size: 75%;"
1518
| colspan="1" | '''Figure 18:''' Detail of the deformed shape with finite element discretization and  distribution of effective plastic strain 
1519
|}
1520
1521
'''Remark'''. The examples presented show that the FIC/FEM formulation is an effective procedure for solving bulk metal forming problems involving full or quasi-incompressible situations. The key advantage of the FIC approach versus more standard mixed FEM formulations is that it provides a natural theoretical framework for equal order finite element interpolations for the velocity and pressure variables, both in the context of implicit and explicit solution schemes. We note the simplicity and effectiveness of the full explicit algorithm as demonstrated in the examples presented. The FIC formulation reproduces also the best feature of the so called stabilized FEM method for incompressible problems, such as the CBS scheme <span id='citeF-20'></span>[[#cite-20|[20,21,25,26,58]]], the pressure gradient operator method <span id='citeF-22'></span>[[#cite-22|[22]]] and the subgrid scale method <span id='citeF-23'></span>[[#cite-23|[23]]] among others.
1522
1523
==8 LAGRANGIAN FLOWS. THE PARTICLE FINITE ELEMENT  METHOD==
1524
1525
===8.1 The Particle Finite Element Method (PFEM)===
1526
1527
The  Lagrangian formulation is  an excellent procedure for treating bulk forming processes involving the interaction of fluids and solids using a unified formulation. An important advantage of the Lagrangian formulation is that both the  motions of the solid  and the fluid  are defined in the same frame of reference and modelled with the some governing equations.
1528
1529
The Lagrangian fluid flow equations can be simply obtained by noting that the velocity of the mesh nodes and that of the fluid particles are the same. Hence the the convective terms vanish in the momentum equations, while the rest of the fluid flow equations remain unchanged. The resulting governing equations have an identical form as those of the Stokes flow problem, with the motion of the flow particles being referred now to a Lagrangian coordinate frame.
1530
1531
The FEM algorithms for solving the Lagrangian flow equations are very similar to those presented for incompressible solids in a previous section. Here we present a particular class of Lagrangian formulation  called the ''particle finite element method'' (PFEM) <span id='citeF-38'></span>[[#cite-38|[38-41]]]. The PFEM treats the mesh nodes in the fluid and solid domains as  dimensionless particles which can freely move an even separate from the main fluid domain representing, for instance, the effect of fluid  drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion.
1532
1533
The quality of the numerical solution   depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.
1534
1535
A typical solution with the PFEM involves the following steps.
1536
1537
<ol>
1538
1539
<li>Identify the external boundaries for both the fluid and solid domains. This is  an essential step as some boundaries (such as the free surface in the fluid) may be severely distorted during the solution process including separation and re-entring of nodes. The Alpha Shape method  is used for the boundary definition (see Section 8.4).  </li>
1540
1541
<li>Discretize the fluid and solid domains with a finite element mesh. For the mesh generation process we use and extended Delaunay technique  [52].  </li>
1542
1543
<li>Solve iteratively the coupled Lagrangian equations of motion for  the fluid and the solid domains. Compute the relevant state variables in both domains at each time step: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid.  </li>
1544
1545
The iterative solution scheme chosen is a particularization of the fractional step algorithm of Section 3.1 for <math display="inline">\alpha =1</math>. In summary the solution steps are the following.
1546
1547
''Step 3.1.'' Compute the predicted velocities (viz, Eq.(31))
1548
1549
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1550
|-
1551
| 
1552
{| style="text-align: left; margin:auto;width: 100%;" 
1553
|-
1554
| style="text-align: center;" | <math>
1555
\tilde{\bf v}^{n+1,i+1} = \bar{\bf v}^n - \Delta t {\bf M}^{-1}_d [{\bf K} \bar{\bf v}^{n} - {\bf G}\bar {\bf p}^{n} -{\bf f}^{n+1}] </math>
1556
|}
1557
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
1558
|}
1559
1560
''Step 3.2.'' Compute <math display="inline">\bar {\textbf p}^{n+1,i+1}</math> from Eq.(32) for <math display="inline">\alpha =1</math>.
1561
1562
''Step 3.3.'' Compute explicitely <math display="inline"> \bar{\textbf v}^{n+1,i+1}</math> from Eq.(33) with <math display="inline">\alpha =1</math>.
1563
1564
''Step 3.4.'' Compute <math display="inline"> \bar{\boldsymbol \pi }^{n+1,i+1}</math> explicitely from Eq.(35).
1565
1566
''Step 3.5.'' Solve for the motion of the solid. This can be performed by integrating the dynamic equations of motion in the solid. Here the algorithm of Section 6 can be used, as it applies to both “compressible” and incompressible materials.
1567
1568
''Step 3.6.'' Estimate a new position of the mesh nodes in terms of the time increment size as
1569
1570
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1571
|-
1572
| 
1573
{| style="text-align: left; margin:auto;width: 100%;" 
1574
|-
1575
| style="text-align: center;" | <math>
1576
{\bf x}_j^{n+1,i+1} = {\bf x}_j^{n}+\bar {\bf u}_j^{n+1,i+1} \Delta t</math>
1577
|}
1578
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
1579
|}
1580
1581
It is important to note that all matrices in Steps 3.1&#8211;3.5 are evaluated at the last predicted configuration <math display="inline">\Omega ^{n+1,i}</math>.
1582
1583
In steps 3.1&#8211;3.6 superindex <math display="inline">i</math> denotes the iteration within each time step.
1584
1585
''Step 3.7.'' Check the convergence of the velocity and pressure fields in the fluid and the displacements, strains and stresses in the solid. If convergence is achieved froze the final position of the mesh nodes and move to the next time increment (step 4), otherwise return to step 3.1 for the next iteration.
1586
1587
<li>Go back to step 1 and repeat the solution process for the next time step. </li>
1588
1589
</ol>
1590
1591
Above algorithm can be found to be analogous to the standard ''updated lagrangian'' scheme typically used in non linear solid mechanics problems <span id='citeF-8'></span>[[#cite-8|[8,57]]].
1592
1593
Despite the motion of the nodes within the iterative process, in general there is no need to regenerate the mesh at each iteration.  ''In the examples presented in the paper the mesh in the fluid domain has been regenerated at each time step''. A cheaper alternative is to generate a new mesh only after a prescribed number of converged time steps, or when  the nodal displacements induce significant geometrical distorsions in some elements.
1594
1595
The boundary conditions are applied as described in Section 3.1.
1596
1597
In the examples presented in the paper the time increment size has been chosen as
1598
1599
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1600
|-
1601
| 
1602
{| style="text-align: left; margin:auto;width: 100%;" 
1603
|-
1604
| style="text-align: center;" | <math>\Delta t =\min (\Delta t_i ) \quad \hbox{with}\quad \Delta t_i ={\vert {\bf v}\vert \over h_i^{\min }} </math>
1605
|}
1606
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
1607
|}
1608
1609
where <math display="inline">h_i^{\min }</math> is the distance between node <math display="inline">i</math> and the closest node in the mesh.
1610
1611
===8.2 Treatment of contact between fluid and solid interfaces===
1612
1613
The  condition of prescribed velocities or pressures at the solid boundaries in the PFEM are  applied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. In some problems it is useful to define a layer of nodes adjacent to the external boundary in the fluid where the condition of prescribed velocity is imposed. These nodes typically remain fixed during the solution process. Contact between liquid particles and the solid boundaries is accounted for by the incompressibility condition which ''naturally prevents the liquid nodes to penetrate into the solid boundaries'' <span id='citeF-38'></span>[[#cite-38|[38-41]]].
1614
1615
===8.3 Generation of a new mesh===
1616
1617
One of the key points for the success of the PFEM  is the fast regeneration of a finite element mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is typically generated using the so called extended Delaunay tesselation (EDT) <span id='citeF-38'></span>[[#cite-38|[38,52]]]. The EDT allows one to generate non standard meshes combining elements of  arbitrary polyhedrical shapes  (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order <math display="inline">n</math>, where <math display="inline">n</math> is the total number of nodes in the mesh. The <math display="inline">C^\circ </math> continuous shape functions of the elements are obtained using the so called meshless finite element interpolation (MFEM) <span id='citeF-53'></span>[[#cite-53|[53]]].
1618
1619
===8.4 Identification of boundary surfaces===
1620
1621
The PFEM requires the correct definition of the boundary domain. Sometimes, boundary nodes are clearly distinguished from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.
1622
1623
Considering that the nodes follow a variable <math display="inline">h(x)</math>  distribution, where <math display="inline">h(x)</math> is the minimum distance between two nodes, the following criterion has been used. For each two nodes (three nodes in 3D) a circle (a sphere in 3D) of radius equal to <math display="inline">\alpha h</math> containing the nodes is plotted. ''All nodes laying on a circle (sphere) with a radius  greater than <math>\alpha h</math>, not containing any other node in the interior are considered as boundary nodes''. In practice, <math display="inline">\alpha </math>  is a parameter close to, but greater than one. This criterion coincides with the Alpha Shape concept <span id='citeF-39'></span>[[#cite-39|[39,55]]].
1624
1625
The method  also allows to identify isolated fluid particles (nodes) outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value. We note that the “flying particles” are in fact “points” ''with no mass'' which motion is followed by integrating the dynamic equations of continuum mechanics for known values of the mass, the initial velocity and acceleration, gravity body forces and a prescribed zero (atmospheric) pressure.
1626
1627
===8.5 Temperature coupled effects===
1628
1629
Thermal-mechanical problems can be effectively treated with the PFEM. This requires solving for the temperature field at each time step, accounting for the couplings  induced by the mechanical equations. The effect of temperature in the mechanical problem is introduced via the constitutive equation in the usual manner.
1630
1631
The form of the heat transfer equation is identical to that of Eq.(70). Recall  that the  Lagrangian formulation eliminates the convective term in the heat transfer equation and hence the FIC stabilization is here not needed <span id='citeF-41'></span>[[#cite-41|[41]]].
1632
1633
===Remark===
1634
1635
The key advantage of the PFEM versus conventional  particle method is that it retains the best feature of the FEM to solve problems in continuum mechanics using a variational approach. A standard finite element mesh is used to discretize in space the problem variables and for solving the governing equations precisely as in the standard FEM. The combination of the lagrangian FEM with the Alpha-Shape approach for identification of the  domain boundary at each time step and the frequent regeneration of the finite element mesh are the distinct features of the PFEM.
1636
1637
==9 APPLICATIONS OF THE PFEM TO METAL FORMING PROCESSES==
1638
1639
===9.1 Sloshing problem===
1640
1641
The sloshing example shown in Figure 19 ilustrates the ability of the PFEM to simulate the flow of liquids within closed cavities with large motions of the free surface. This feature of the PFEM is essential to model the filling of moulds as shown in the next examples. More applications of the  PFEM to sloshing problemas are reported in <span id='citeF-39'></span>[[#cite-39|[39]]].
1642
1643
<div id='img-19'></div>
1644
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1645
|-
1646
|[[File:Draft_Samper_805805629_3102_Fig19.png|501px|PFEM results for a large amplitude sloshing problem]]
1647
|- style="text-align: center; font-size: 75%;"
1648
| colspan="1" | '''Figure 19:''' PFEM results for a large amplitude sloshing problem
1649
|}
1650
1651
===9.2 Filling of moulds===
1652
1653
Figure 20 shows the results of the filling of a 2D mould cavity with a powder material using the PFEM. Note that the essential feature of the filling process are well reproduced. The finite element mesh used for the computation in a certain instant is shown in Figure 20b. This ilustrates the fact that the PFEM is, in fact, a blending of particle and finite element procedure. Other applications of the PFEM to powder compaction problems are reported in <span id='citeF-57'></span>[[#cite-57|[57]]].
1654
1655
Figure 21 shows the different stages of the 3D simulation of the filling of a mould with a casted metal using the PFEM. The pictures show the evolution of the casting during the filling process. The metal flows through the vertical mould first and then fills the lower and upper lateral ducts, as expected. The gray colour intensities denote the different temperatures in the melt with a whiler colour indicating the cooler region. Note that the metal cools down  as it progresses within the mould.
1656
1657
The example in Figure 22 shows the simulation of the casting of a mechanical part using the PFEM. Note that the mesh discretizing the casted region is progressively generated as the mould is filled. This is a distinct feature of the PFEM which distinguishes this method from more standard VOF techniques using Eulerian and ALE approaches.
1658
1659
<div id='img-20'></div>
1660
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1661
|-
1662
|[[Image:Draft_Samper_805805629-Figure19_1.png|390px|]]
1663
|[[Image:Draft_Samper_805805629-Figure19_2.png|390px|]]
1664
|-
1665
|[[Image:Draft_Samper_805805629-Figure19_3.png|390px|]]
1666
|[[Image:Draft_Samper_805805629-Figure19_4.png|390px|]]
1667
|-
1668
|[[Image:Draft_Samper_805805629-Figure19_5.png|390px|]]
1669
|[[Image:Draft_Samper_805805629-Figure19_6.png|390px|]]
1670
|-
1671
|colspan="2" | (a)
1672
|-
1673
|colspan="2" | [[Image:Draft_Samper_805805629-cuadro25mesh.png|470px|]]
1674
|-
1675
|colspan="2" | (b)
1676
|- style="text-align: center; font-size: 75%;"
1677
| colspan="2" | '''Figure 20:''' (a) Filling of a 2D mould with a powder material using the PFEM. (b) Finite element mesh used for the computations in a certain instant. The circles indicate the nodes at the external and internal free surfaces
1678
|}
1679
1680
1681
<div id='img-22'></div>
1682
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1683
|-
1684
|[[Image:Draft_Samper_805805629-fig1.png|300px|]]
1685
|[[Image:Draft_Samper_805805629-fig2.png|300px|]]
1686
|-
1687
|[[Image:Draft_Samper_805805629-fig3.png|300px|]]
1688
|[[Image:Draft_Samper_805805629-fig4.png|300px|Filling of a 3D mould with a casted metal using the PFEM. Colours indicate the temperature values]]
1689
|- style="text-align: center; font-size: 75%;"
1690
| colspan="2" | '''Figure 21:''' Filling of a 3D mould with a casted metal using the PFEM. Colours indicate the temperature values
1691
|}
1692
1693
<div id='img-23'></div>
1694
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1695
|-
1696
|[[Image:Draft_Samper_805805629-picture1.png|240px|]]
1697
|[[Image:Draft_Samper_805805629-picture50.png|240px|]]
1698
|[[Image:Draft_Samper_805805629-picture150.png|240px|]]
1699
|-
1700
|[[Image:Draft_Samper_805805629-picture300.png|240px|]]
1701
|[[Image:Draft_Samper_805805629-picture600.png|240px|]]
1702
|[[Image:Draft_Samper_805805629-picture801.png|240px|Simulation of a casting process using the PFEM]]
1703
|- style="text-align: center; font-size: 75%;"
1704
| colspan="3" | '''Figure 22:''' Simulation of a casting process using the PFEM
1705
|}
1706
1707
<div id='img-24'></div>
1708
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1709
|-
1710
|[[Image:Draft_Samper_805805629-medit.500.png|351px|(cont.) Cut through the mesh]]
1711
|- style="text-align: center; font-size: 75%;"
1712
| colspan="1" | '''Figure 22:''' (cont.) Cut through the mesh
1713
|}
1714
1715
===9.3 Mixing problems===
1716
1717
The PFEM is particularly suited for analyzing problems involving mixing of two fluids or that of a collection of particles in a fluid. This is very useful to model a wide class of material forming processes <span id='citeF-5'></span>[[#cite-5|[5,7]]].
1718
1719
A first example of this kind is shown in Figure 23. The problem represents the mixing of two fluids of different density. The PFEM allows to track the position of the fluid particles of lower density which mix within the second fluid.
1720
1721
<div id='img-25'></div>
1722
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1723
|-
1724
|[[Image:Draft_Samper_805805629-mixing1.png|600px|]]
1725
|-
1726
|[[Image:Draft_Samper_805805629-mixing2.png|600px|]]
1727
|-
1728
| [[Image:Draft_Samper_805805629-mixing3.png|350px|Analysis of the mixing of two fluids of different density with the PFEM]]
1729
|- style="text-align: center; font-size: 75%;"
1730
| '''Figure 23:''' Analysis of the mixing of two fluids of different density with the PFEM
1731
|}
1732
1733
Figure 24 shows a last example of application of the PFEM to the mixing of a collection of particles within a container filled with a fluid of a higher density. Initially the particles are thrown into the container and mix within the fluid as shown. As time evolves the particles move up naturally towards the surface of the fluid due to their lower density.
1734
1735
The last two  examples clearly show the possibilities of the PFEM for analysis of  material mixing situations. Extensions of the approach accounting for thermal-coupled effects and chemical reactive situations are straightforward <span id='citeF-41'></span>[[#cite-41|[41]]].
1736
1737
<div id='img-26'></div>
1738
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1739
|-
1740
|[[Image:Draft_Samper_805805629-Fig_mixing.png|578px|Mixing of particles in a fluid. Evolution of the particles during the mixing process]]
1741
|- style="text-align: center; font-size: 75%;"
1742
| colspan="1" | '''Figure 24:''' Mixing of particles in a fluid. Evolution of the particles during the mixing process
1743
|}
1744
1745
==10 CONCLUDING REMARKS==
1746
1747
The finite calculus form of the balance equations  in mechanics is a good starting point for deriving stabilized finite element  methods for solving a variety of bulk metal forming problems involving fluid flow  and large deformation of solids using Euler and  Lagrangian descriptions. Free surface  effects and coupled thermal situations  can be accounted for in a straightforward manner within the general   algorithm.  A fractional step scheme with intrinsic stabilization properties has been described. Applications of the Eulerian and Lagrangian methods presented range from mould filling, solidification and cooling situations in casting processes to forging and extrussion problems, among many others. The Lagrangian PFEM formulation is very effective for the analysis of bulk metal forming problems and material design processes involving filling of moulds,  material mixing situations and complex fluid-structure interactions.
1748
1749
==ACKNOWLEDGEMENT==
1750
1751
The authors thank the company QUANTECH ATZ, SA (www.quantech.es) for providing the code VULCAN for solution of the casting problem presented. Thanks are also given to MIKALOR SA for providing data for the hose clamp example.
1752
1753
==BIBLIOGRAPHY==
1754
1755
<div id="cite-1"></div>
1756
'''[[#citeF-1|[1]]]'''  J.L. Chenot and E. Oñate, ''Modelling of Metal Forming Processes'', Kluwer Academic, Dordrecht, 1988.
1757
1758
<div id="cite-2"></div>
1759
'''[[#citeF-2|[2]]]'''  J.F.T. Pittman, O.C. Zienkiewicz, R.D. Wood and J.M. Alexander (Eds), '' Numerical Analysis of Forming Processes'', Wiley, Chichester, 1984.
1760
1761
<div id="cite-3"></div>
1762
'''[[#citeF-3|[3]]]'''  J. Huetnik and  F.P.T. Baaijens (Eds), ''Simulation of Materials Processing: Theory,  Methods and Applications, NUMIFORM 98'',  Enschede, Países Bajos,  A.A. Balkema, Rotterdam, 22&#8211;5 Jun. 1998.
1763
1764
<div id="cite-4"></div>
1765
'''[[#citeF-4|[4]]]'''  M. Pietrzyk, J. Kusiak, J. Majta, P. Hartley and I. Pillinger, ''Metal Forming 2000'', A.A. Balkema, Rotterdam, 2000.
1766
1767
<div id="cite-5"></div>
1768
'''[[#citeF-5|[5]]]'''  K. Mori, ''Simulation of Material Processing: Theory, Methods and Applications'', A.A. Balkema, Lisse, 2001.
1769
1770
<div id="cite-6"></div>
1771
'''[[#citeF-6|[6]]]'''  D.R.J. Owen, E. Oñate and B. Suarez (Eds.), ''Computational Plasticity VII. Fundamentals and Applications'', CIMNE, Barcelona 2003.
1772
1773
<div id="cite-7"></div>
1774
'''[[#citeF-7|[7]]]'''  S. Ghosh, J.C. Castro and J.K. Lee, ``Material Processing and Design: Simulation and ApplicatOñate,ions, NUMIFORM 2004, American Institute of Physics, 2004.
1775
1776
<div id="cite-8"></div>
1777
'''[[#citeF-8|[8]]]''' O.C. Zienkiewicz  and  R.C. Taylor,''The finite element method'', 5th Edition, 3 Volumes, Butterworth&#8211;Heinemann, 2000.
1778
1779
<div id="cite-9"></div>
1780
'''[[#citeF-9|[9]]]''' C. Hirsch, ''Numerical computation of internal and external flow'', J. Wiley, Vol. 1 1988, Vol. 2, 1990.
1781
1782
<div id="cite-10"></div>
1783
'''[[#citeF-10|[10]]]''' J. Peraire, K. Morgan  and J. Peiro, “The simulation of 3D incompressible flows using unstructured grids”, In ''Frontiers of Computational Fluid Dynamics'', Caughey DA and Hafez MM. (eds), Chapter 16, J. Wiley, 1994.
1784
1785
<div id="cite-11"></div>
1786
'''[[#citeF-11|[11]]]''' M. Storti, N. Nigro  and S.R. Idelsohn, “Steady state incompressible flows using explicit schemes with an optimal local preconditioning”, ''Computer Methods in Applied Mechanics and Engineering'', '''124''', 231&#8211;252, 1995.
1787
1788
<div id="cite-12"></div>
1789
'''[[#citeF-12|[12]]]''' J.C. Heinrich, P.S. Hayakorn and O.C. Zienkiewicz, “An upwind finite element scheme for two dimensional convective transport equations”, ''Int. J. Num. Meth. Engng.'', '''11''', 131&#8211;143, 1977.
1790
1791
<div id="cite-13"></div>
1792
'''[[#citeF-13|[13]]]''' A. Brooks and T.J.R. Hughes, “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.
1793
1794
<div id="cite-14"></div>
1795
'''[[#citeF-14|[14]]]''' T.J.R. Hughes and M. Mallet, “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''', pp. 305&#8211;328,  1986.
1796
1797
<div id="cite-15"></div>
1798
'''[[#citeF-15|[15]]]'''  M.A. Cruchaga and E. Oñate, “A finite element formulation for incompressible flow problems using a generalized streamline operator”, ''Comput. Methods in Appl. Mech. Engrg.'', '''143''', 49&#8211;67, 1997.
1799
1800
<div id="cite-16"></div>
1801
'''[[#citeF-16|[16]]]'''   M.A. Cruchaga and E. Oñate, “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.
1802
1803
<div id="cite-17"></div>
1804
'''[[#citeF-17|[17]]]''' T.J.R. Hughes, L.P. Franca and G.M. Hulbert, “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''',  pp. 173&#8211;189, 1989.
1805
1806
<div id="cite-18"></div>
1807
'''[[#citeF-18|[18]]]''' J. Donea, “A Taylor-Galerkin method for convective transport problems”, ''Int. J. Num. Meth. Engng.'', '''20''', 101&#8211;119, 1984.
1808
1809
<div id="cite-19"></div>
1810
'''[[#citeF-19|[19]]]''' R. Löhner, K. Morgan and O.C. Zienkiewicz, “The solution of non-linear hyperbolic equation systems by the finite element method”, ''Int. J. Num. Meth. in Fluids'', '''4''', 1043, 1984.
1811
1812
<div id="cite-20"></div>
1813
'''[[#citeF-20|[20]]]'''  R. Codina, M. Vazquez and O.C. Zienkiewicz,  “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.
1814
1815
<div id="cite-21"></div>
1816
'''[[#citeF-21|[21]]]'''  R. Codina and O.C. Zienkiewicz, “CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter”, ''Communications in Numerical Methods in Engineering'', '''18 (2)''', 99&#8211;112, 2002.
1817
1818
<div id="cite-22"></div>
1819
'''[[#citeF-22|[22]]]'''  R. Codina and O. Soto, “A numerical model to track two-fluid interfaces based on a stabilized finite element method and the level set technique”, ''Int. J. Num. Meth. in Fluids'', '''4''', 293&#8211;301, (2002).
1820
1821
<div id="cite-23"></div>
1822
'''[[#citeF-23|[23]]]'''  R. Codina, “Stabilized finite element approximation of transient incompressible flows using orthogonal subscales”, ''Comput. Methods Appl. Mech. Engrg.'', '''191''', 4295&#8211;4321, 2002.
1823
1824
<div id="cite-24"></div>
1825
'''[[#citeF-24|[24]]]'''  J. Donea and A. Huerta, ''"Finite element method for flow problems"'', J. Wiley, 2003.
1826
1827
<div id="cite-25"></div>
1828
'''[[#citeF-25|[25]]]'''  O.C. Zienkiewicz, J. Rojek, R.L. Taylor and M. Pastor, “Triangles and tetrahedra in explicity dynamic codes for solids”, ''Int. J. Num. Meth. Engng.'', '''43''', 565&#8211;583, 1998.
1829
1830
<div id="cite-26"></div>
1831
'''[[#citeF-26|[26]]]'''  J. Rojek, O.C. Zienkiewicz,  E. Oñate and E. Postek, “Advances in FE explicit formulations for simulation of metal forming processes”, ''J. of Material Processing Technol.'', '''119''', 41&#8211;47, 2001.
1832
1833
<div id="cite-27"></div>
1834
'''[[#citeF-27|[27]]]'''   M. Chiumenti, Q. Valverde, C. Agelet de Saracibar and M. Cervera, “A stabilized formulation for incompressible elasticity using linear displacement and pressure interpolations”, ''Computer Methods in Applied Mechanics and Engineering'', '''191''', ('''46'''), 5253&#8211;5264, (2002).
1835
1836
<div id="cite-28"></div>
1837
'''[[#citeF-28|[28]]]'''  M. Cervera, M. Chiumenti, Q. Valverde and C. Agelet de Saracibar, “Mixed linear/linear simplicial elements for incompressible elasticity and plasticity”, ''Comput. Meth. Appl. Mech. Engrg.'', '''192''' ('''49-50'''), 5249&#8211;5263, (2003).
1838
1839
<div id="cite-29"></div>
1840
'''[[#citeF-29|[29]]]'''  M. Chiumenti, Q. Valverde, C. Agelet de Saracibar  and M. Cervera, “A stabilized formulation for incompressible plasticity using linear triangles and tetrahedra”, ''Int. Journal of Plasticity'', '''20''' ('''8-9'''), 1487&#8211;1504, (2004).
1841
1842
<div id="cite-30"></div>
1843
'''[[#citeF-30|[30]]]''' E. Oñate, Derivation of stabilized equations for advective-diffusive transport and fluid flow problems, ''Comput. Meth. Appl. Mech. Engng.'', Vol. '''151''', pp. 233&#8211;267, (1998).
1844
1845
<div id="cite-31"></div>
1846
'''[[#citeF-31|[31]]]'''  E. Oñate, “Possibilities of finite calculus in computational mechanics”,  ''Int. J. Num. Meth. Engng.'', '''60'''(1), 255&#8211;281, 2004.
1847
1848
<div id="cite-32"></div>
1849
'''[[#citeF-32|[32]]]'''  E. Oñate, “A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation”, ''Comp. Meth. Appl. Mech. Engng.'', '''182''', 1&#8211;2, 355&#8211;370, 2000.
1850
1851
<div id="cite-33"></div>
1852
'''[[#citeF-33|[33]]]'''  E. Oñate and J. García, “A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation”, ''Comput. Meth. Appl. Mech. Engrg.'', '''191''', 635&#8211;660, (2001).
1853
1854
<div id="cite-34"></div>
1855
'''[[#citeF-34|[34]]]'''  E. Oñate, J. García, S.R. Idelsohn and F. Del Pin. FIC formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches, Submitted to ''Comput. Methods Appl. Mech. Engng.'', July 2004.
1856
1857
<div id="cite-35"></div>
1858
'''[[#citeF-35|[35]]]'''  E. Oñate, J. Arteaga, J, García and R. Flores, “Error estimation and mesh adaptivity in incompressible viscous flows using a residual power approach”, Submitted to ''Comput. Meth. Appl. Mech. Engrg.'', 2004.
1859
1860
<div id="cite-36"></div>
1861
'''[[#citeF-36|[36]]]'''  E. Oñate, R.L. Taylor, O.C. Zienkiewicz and J. Rojek, “A residual correction method based on finite calculus”, ''Engineering Computations'', '''20''', (5/6), 629&#8211;658, 2003.
1862
1863
<div id="cite-37"></div>
1864
'''[[#citeF-37|[37]]]'''   E. Oñate, J. Rojek, R.L. Taylor and O.C. Zienkiewicz, “Finite calculus formulation for incompressible solids using linear triangles and tetrahedra”, ''Int. J. Num. Meth. Engng.'',  '''59''', 1473&#8211;1500, 2004.
1865
1866
<div id="cite-38"></div>
1867
'''[[#citeF-38|[38]]]'''  S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, ''Computer and Structures'', '''81''', 655&#8211;671, 2003.
1868
1869
<div id="cite-39"></div>
1870
'''[[#citeF-39|[39]]]'''  E. Oñate, S.R. Idelsohn, F. Del Pin and R. Aubry,  “The Particle finite element method. An overview”, ''Int. J. Comput. Meth.'', '''1''' ('''2'''), 267&#8211;307, 2004.
1871
1872
<div id="cite-40"></div>
1873
'''[[#citeF-40|[40]]]'''  S.R. Idelsohn, E. Oñate and F. Del Pin, “The Particle Finite Element Method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves”, ''Int. J. Num. Meth. Engng.'', (Submitted), 2004.
1874
1875
<div id="cite-41"></div>
1876
'''[[#citeF-41|[41]]]'''  R. Aubry, S.R. Idelsohn  and E. Oñate, “Particle finite element method in fluid mechanics including thermal convection-diffusion”. ''Computer & Structures'', submitted, 2004.
1877
1878
<div id="cite-42"></div>
1879
'''[[#citeF-42|[42]]]'''  E. Oñate, C. Sacco and  S. Idelsohn, “A finite point method for incompressible flow problems”, ''Computing and Visualization in Science'', '''2''', 67&#8211;75, 2000.
1880
1881
<div id="cite-43"></div>
1882
'''[[#citeF-43|[43]]]'''  O.C. Zienkiewicz, P.C. Jain y E. Oñate, Flow of solids during forming and extrusion: some aspects of numerical solutions, ''Int. J. Solids Struct.'', '''14''', 15&#8211;38, 1978.
1883
1884
<div id="cite-44"></div>
1885
'''[[#citeF-44|[44]]]'''  O.C. Zienkiewicz, E. Oñate and J.C. Heinrich,  A general formulation for coupled thermal flow of metals using finite elements, ''Int. J. Num. Meth. Eng.'', ''' 17''', 1497&#8211;514, 1981.
1886
1887
<div id="cite-45"></div>
1888
'''[[#citeF-45|[45]]]'''  R. Codina, “Pressure stability in fractional step finite element method for incompressible flows”, ''J. Comput. Physics'', '''170''', 112&#8211;140, 2001.
1889
1890
<div id="cite-46"></div>
1891
'''[[#citeF-46|[46]]]''' E. Oñate and M. Manzan, 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.
1892
1893
<div id="cite-47"></div>
1894
'''[[#citeF-47|[47]]]'''  Codina, R., Schaefer, U., and Oñate, E., “Mould filling simulation using finite elements”, ''Int. J. Num. Meth. Engng.'', '''4''', 291&#8211;310, (1994).
1895
1896
<div id="cite-48"></div>
1897
'''[[#citeF-48|[48]]]''' Codina, R., and Soto, O., “Finite element simulation of the filling of thin moulds”, ''Int. J. Num. Meth. Engng.'', '''46''', 1559&#8211;1573, (1999).
1898
1899
<div id="cite-49"></div>
1900
'''[[#citeF-49|[49]]]'''  C. Agelet de Saracibar, M. Cervera and M. Chiumenti, “ On the Formulation of coupled thermoplastic problems with phase-change”, ''International Journal of Plasticity'', Vol. '''15''', 1&#8211;34, (1999).
1901
1902
<div id="cite-50"></div>
1903
'''[[#citeF-50|[50]]]'''  M. Cervera, C. Agelet de Saracibar and M. Chiumenti, “ Thermo-mechanical analysis of industrial solidification processes”, ''Int. Journal for Num. Methods in Engineering'', Vol. '''46''', 1575-1591, (1999).
1904
1905
<div id="cite-51"></div>
1906
'''[[#citeF-51|[51]]]'''  C. Agelet de Saracibar, M. Cervera and M. Chiumenti, “On the constitutive modeling of coupled thermo-mechanical phase-change problems”, ''International Journal of Plasticity'', Vol. '''17''', 1565-1622, (2001).
1907
1908
<div id="cite-52"></div>
1909
'''[[#citeF-52|[52]]]'''  N. Calvo, S.R. Idelsohn and E. Oñate, “The extended Delaunay tesselation”, ''Engineering Computations'', '''20''' (5/6), 583&#8211;600, 2003.
1910
1911
<div id="cite-53"></div>
1912
'''[[#citeF-53|[53]]]''' S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”,  ''Int. J. Num. Meth. Engng.'', Vol. '''58''', pp. 893- 912, 2003.
1913
1914
<div id="cite-54"></div>
1915
'''[[#citeF-54|[54]]]'''  MARC/Autoforge User's Guide 2.2. 1999.
1916
1917
<div id="cite-55"></div>
1918
'''[[#citeF-55|[55]]]'''  H. Edelsbrunner  and E.P. Mucke, E.P., “Three dimensional alpha shapes”, ''ACM Trans. Graphics'', '''13''', 43&#8211;72, 1999.
1919
1920
<div id="cite-56"></div>
1921
'''[[#citeF-56|[56]]]'''  VULCAN. Finite element system for casting simulation, Quantech ATZ, Barcelona, Spain, (2003), (www.quantech.es).
1922
1923
<div id="cite-57"></div>
1924
'''[[#citeF-57|[57]]]'''  J.C. Cante, J. Oliver and C. González, “Powder forming simulation with the particle finite element method”, ''Computational Mechanics'', Z.H. Yao, M.W. Yuan and W.X. Zhong (Eds.), CD-Rom Proceedings of the WCCM VI, Beijing, September 5&#8211;10, 2004.
1925
1926
<div id="cite-58"></div>
1927
'''[[#citeF-58|[58]]]''' J. Rojek, E. Oñate and R.L. Taylor, “CBS based stabilization in explicit solid mechanics”. ''Int. J. Num. Mech. Engng.'' Accepted for publication. October 2005.
1928
1929
==APPENDIX==
1930
1931
In Eqs.(24)
1932
1933
<span id="eq-A.1"></span>
1934
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1935
|-
1936
| 
1937
{| style="text-align: left; margin:auto;width: 100%;" 
1938
|-
1939
| style="text-align: center;" | <math>{\bf H}={\bf A}+{\bf K}+\hat {\bf K}  </math>
1940
|}
1941
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
1942
|}
1943
1944
If we denote the node indexes with superscripts <math display="inline">a,b</math>, the space indexes with subscripts <math display="inline">i,j</math>, the element contributions to the components of the arrays involved in Eqs.(24) and (A.1) are
1945
1946
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1947
|-
1948
| 
1949
{| style="text-align: left; margin:auto;width: 100%;" 
1950
|-
1951
| style="text-align: center;" | <math>\begin{array}{l} \displaystyle M_{ij}^{ab}= \left(\int _{\Omega ^e} \rho N^a N^b d\Omega \right)\delta _{ij} \quad ,\quad A_{ij}^{ab}= \left(\int _{\Omega ^e} \rho N^a ({\textbf u}^T {\nabla } N^b) d\Omega \right)\delta _{ij}\\ \\ 
1952
\displaystyle {K}_{ij}^{ab} = \left(\int _{\Omega ^e} \mu {\boldsymbol \nabla }^T N^a{\boldsymbol \nabla }N^b  d\Omega \right)\delta _{ij} \quad ,\quad {\boldsymbol \nabla } = \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2},{\partial  \over \partial x_3}\right]^T\\ \\ 
1953
\displaystyle \hat{K}_{ij}^{ab} = \left({1\over 2} \int _{\Omega ^e} ({\textbf h}^T {\boldsymbol \nabla } N^a)(\rho {\textbf u}^T {\boldsymbol \nabla } N^b)d\Omega \right)\delta _{ij}\quad ,\quad {G}_{i}^{ab}=\int _{\Omega ^e} {\partial N^a \over \partial x_i}N^b d\Omega \\ \\ 
1954
\displaystyle {\textbf C}= \left[\begin{matrix}{\textbf C}_1\\ {\textbf C}_2\\ {\textbf C}_3\\\end{matrix}\right]\quad ,\quad {C}_{ij}^{ab}=\left({1\over 2} \int _{\Omega ^e} [{\textbf h}^T {\boldsymbol \nabla } N^a]N^bd\Omega \right)\delta _{ij}\\ \\ 
1955
\displaystyle \hat L^{ab}= \int _{\Omega ^e} ({\boldsymbol \nabla }^T N^a) [\tau ] {\boldsymbol \nabla } N^b d\Omega \quad ,\quad [\tau ]= \left[\begin{matrix}\tau _1 &0 &0 \\ 0 & \tau _2&0\\ 0&0& \tau _3\\\end{matrix}\right]\\ \\ 
1956
\displaystyle {\textbf Q}= [{\textbf Q}_1,{\textbf Q}_2,{\textbf Q}_3] \quad ,\quad  \displaystyle Q_{i}^{ab} = \int _{\Omega ^e}\tau _i {\partial N^a \over \partial x_i} N^b d\Omega \quad \quad \hbox{no sum in }i\\ \\ 
1957
\displaystyle \hat{\textbf C}= [\hat{\textbf C}_1,\hat{\textbf C}_2,\hat{\textbf C}_3] \quad ,\quad  \displaystyle \hat{C}_1^{ab}=\hat{C}_2^{ab}=\hat{C}_3^{ab} = \int _{\Omega ^e} \rho ^2 N^a ({\textbf u}^T {\boldsymbol \nabla }N^b)d\Omega \\ \\ 
1958
\displaystyle  \hat {M}^{ab}_{ij}= \left(\int _{\Omega ^e} \tau _i N^a N^b d\Omega \right)\delta _{ij}\quad ,\quad  \displaystyle {f}_i^a = \int _{\Omega ^e} N^a f_i d\Omega + \int _{\Gamma ^e} N^a t_i d\Gamma  \end{array}</math>
1959
|}
1960
|}
1961
1962
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1963
|-
1964
| 
1965
{| style="text-align: left; margin:auto;width: 100%;" 
1966
|-
1967
| style="text-align: center;" | 
1968
|}
1969
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
1970
|}
1971
1972
It is understood that all the arrays are matrices (except <math display="inline">{f}</math>, which is a vector) whose components are obtained by grouping together the left indexes in the previous expressions (<math display="inline">a</math> and possibly <math display="inline">i</math>) and the right indexes (<math display="inline">b</math> and possibly <math display="inline">j</math>).
1973
1974
Note that the stabilization matrix <math display="inline">\hat{K}</math> in Eq.([[#eq-A.1|A.1]]) adds additional orthotropic diffusivity terms of value <math display="inline">\rho \displaystyle{h_ku_l\over 2}</math>.
1975

Return to Onate Rojek et al 2006a.

Back to Top

Document information

Published on 01/01/2006

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

Document Score

0

Times cited: 6
Views 31
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?