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
==A particle finite element method for  analysis of industrial forming processes==
2
3
'''E. Oñate''', '''A. Franci''', '''J.M. Carbonell'''
4
5
Centre Internacional de Metodes Numerics a l'Enginyeria - CIMNE, Barcelona, Spain<br />
6
Universitat Politècnica de Catalunya (UPC), Barcelona, Spain<br />
7
RR
8
9
10
11
==Abstract==
12
13
We present a generalized Lagrangian formulation for analysis of industrial forming processes involving thermally coupled interactions between deformable continua. The governing equations for the deformable bodies are written in a unified manner that holds both for fluids and solids. The success of the formulation lays on  a  residual-based  expression of the mass conservation equation obtained using the Finite Calculus (FIC) method that provides the necessary stability for quasi/fully incompressible situations. The  governing equations are discretized with the FEM via a mixed formulation using simplicial elements with equal linear interpolation for the velocities, the pressure and the temperature. The merits of the formulation  are demonstrated in the solution of 2D and 3D thermally-coupled forming processes using the Particle Finite Element Method (PFEM).
14
15
'''keywords''': Updated Lagrangian formulation, finite element method, incompressible fluids, consistent tangent matrix, mixed formulation, stabilized method
16
17
18
==1 INTRODUCTION==
19
20
There is a large number of  forming manufacturing processes in industry that involve the interaction of highly deformable continua, including viscous fluids and solids that undergo large deformations. Typically these problems also include coupled thermal effects between the interacting bodies during the forming process.
21
22
A number of numerical methods has been developed over the years for the simulation of    industrial forming processes. The different methods can be classified as those based on a pure fluid mechanics formulation, i.e. the so-called flow approach developed around  the 1980's <span id='citeF-6'></span><span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-44'></span><span id='citeF-45'></span><span id='citeF-48'></span>[[#cite-6|[6,39,40,44,45,48]]] and those based on a solid mechanics formulation <span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-19'></span><span id='citeF-31'></span><span id='citeF-47'></span>[[#cite-11|[11,12,19,31,47]]].
23
24
In this work we present a generalized continuum mechanics formulation for the analysis of industrial forming processes that involve thermally coupled interactions between  deformable continua (both fluids and solids).
25
26
The motion of the deformable bodies (either a solid or a fluid) is described using a  Lagrangian mixed formulation with the velocities, the pressure and the temperature as the unknown variables. Appropriate constitutive laws for the different constituent materials in the different parts of the analysis domain  are used. The mixed  Lagrangian formulation allows us to model the non linear interaction (including coupled thermal effects) between the different bodies involved in the forming process in a unified manner.
27
28
The continuum mechanics equations are solved with the Particle Finite Element Method (PFEM, www.cimne.com/pfem). The PFEM treats the mesh nodes in the analysis domain as virtual particles that model the behaviour of the underlying  can freely move and even separate from the  domain representing, for instance, the effect of water drops or cutting particles in drilling problems and machining. A mesh connects the nodes discretizing the continuum domain where the governing equations are solved using a stabilized FEM. Examples of application of PFEM to problems in fluid and solid mechanics including fluid-structure interaction (FSI) situations and coupled thermal flow problems can be found in <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-9'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-20'></span><span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-34'></span><span id='citeF-35'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-43'></span>[[#cite-1|[1,2,4,5,9,13,14,15,16,17,18,20,26,27,34,35,37,38,43]]]. Early attempts of the using the PFEM for solving metal forming problems were reported in <span id='citeF-20'></span><span id='citeF-30'></span>[[#cite-20|[20,30]]].
29
30
In Lagrangian analysis procedures (such as  PFEM) the motion of the virtual particles (i.e. the nodes) and of physical particles is tracked during the transient solution. Hence, the convective terms vanish in the momentum and heat transfer equations for the continuum part of the domain and no numerical stabilization is needed for treating those terms. Another source of numerical instability however remains in Lagrangian formulation, namely the treatment of incompressibility using equal order FEM interpolation for the velocities and the pressure. In this work a stabilized form of the discretized mass conservation equation is obtaining using the Finite Calculus (FIC) method proposed by Oñate et al. <span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-30'></span><span id='citeF-36'></span>[[#cite-21|[21,22,23,24,25,28,29,30,36]]] that has excellent mass preservation features.
31
32
The generalized Lagrangian approach here described follows the ideas presented in <span id='citeF-15'></span>[[#cite-15|[15]]] where a PFEM formulation for treating fluids and solids in a unified  manner was proposed. These ideas have been extended by Franci, Oñate and Carbonell <span id='citeF-10'></span>[[#cite-10|[10]]] using a stabilized PFEM formulation based on the FIC method.
33
34
The lay-out of the paper is the following. In the next section we present the basic equations for conservation of linear momentum, mass conservation and heat transfer for a  continuum using a generalized Lagrangian framework. Next we derive the stabilized FIC form of the mass conservation equation. An incompressible continuum can be considered as the limit case of the compressible problem. Then the mixed finite element discretization using simplicial element with equal order approximation for the velocity,  the pressure and the temperature is presented and the relevant matrices and vectors of the discretized problem are given.  Details of the implicit solution of the Lagrangian FEM equations in time using an updated Lagrangian approach and a Newton-Raphson type iterative scheme are presented. The relevance of the bulk stiffness terms in the tangent matrix  for enhancing the convergence and accuracy of the iterative solution  for problems involving quasi and fully incompressible continua is discussed. The basic steps of the PFEM are described.
35
36
The efficiency and general applicability of the PFEM  are verified by solving a set of thermally coupled industrial forming problems in two (2D) and three (3D) dimensions. The excellent performance of the numerical method proposed for all the problems analyzed is highlighted.
37
38
<div id='img-1'></div>
39
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
40
|-
41
|[[Image:draft_Samper_718004671-Figura1.png|600px|Initial (t =⁰t), current (t=ⁿt) and updated (t=ⁿ⁺¹ t) configurations of a  domain V with Neumann (Γₜ) and Dirichlet (Γ<sub>v</sub>) boundaries. Trajectory of a material point i and velocity  (v<sub>i</sub>) and displacement (u<sub>i</sub>) vectors of the point at each configuration. ⁿu, ⁿ⁺¹u and ∆u denote schematically the trajectories of the overall domain between two configurations.]]
42
|- style="text-align: center; font-size: 75%;"
43
| colspan="1" | '''Figure 1:''' Initial (<math>t ={}^0 t</math>), current (<math>t={}^n t</math>) and updated (<math>t={}^{n+1} t</math>) configurations of a  domain <math>V</math> with Neumann (<math>\Gamma _t</math>) and Dirichlet (<math>\Gamma _v</math>) boundaries. Trajectory of a material point <math>i</math> and velocity  (<math>{v}_i</math>) and displacement (<math>{u}_i</math>) vectors of the point at each configuration. <math>{}^n {u}</math>, <math>{}^{n+1}{u}</math> and <math>\Delta {u}</math> denote schematically the trajectories of the overall domain between two configurations.
44
|}
45
46
==2 GOVERNING EQUATIONS FOR A LAGRANGIAN CONTINUUM==
47
48
We will consider a domain containing a deformable material (either a fluid or a solid) which evolves in time due to external and internal forces and prescribed velocities and thermal conditions from an ''initial configuration'' at time <math display="inline">t={}^0 t</math> (typically <math display="inline">{}^0t =0</math>) to a ''current configuration'' at time <math display="inline">t={}^n t</math>. The volume <math display="inline">V</math> and its boundary <math display="inline">\Gamma </math> at the initial and current configurations are denoted as (<math display="inline">{}^0 V, {}^0 \Gamma </math>) and (<math display="inline">{}^n V, {}^n \Gamma </math>), respectively. The goal is to find the domain that the material  occupies, as well as the velocities, strain rates and stresses (the deviatoric stresses and the pressure) and the temperature in the so-called ''updated configuration''  at time <math display="inline">{}^{n+1} t= {}^n t+\Delta t</math> (Figure 1). In the following a left superindex denotes the configuration where the variable is computed.
49
50
===2.1 Momentum equations===
51
52
The equations of conservation of linear momentum for a deformable continuum are written in the Lagrangian description as <span id='citeF-3'></span>[[#cite-3|[3]]]
53
54
<span id="eq-1"></span>
55
{| class="formulaSCP" style="width: 100%; text-align: left;" 
56
|-
57
| 
58
{| style="text-align: left; margin:auto;" 
59
|-
60
| style="text-align: center;" | <math>\rho \frac{Dv_i}{Dt}-{\partial {}^{n+1}\sigma _{ij} \over \partial {}^{n+1}x_j} - {}^{n+1} b_i=0\quad , \quad i,j=1,\cdots , n_s \quad  \hbox{in } {}^{n+1}V  </math>
61
|}
62
| style="width: 5px;text-align: right;" | (1)
63
|}
64
65
In Eq.([[#eq-1|1]]), <math display="inline">{}^{n+1}V</math> is the analysis domain in the ''updated configuration'' at time <math display="inline">^{n+1}t</math> with boundary <math display="inline">{}^{n+1}\Gamma </math>, <math display="inline">v_i</math> and <math display="inline">b_i</math> are the velocity and body force components  along the <math display="inline">i</math>th Cartesian axis, <math display="inline">\rho </math> is the density, <math display="inline">n_s</math> is the number of space dimensions (i.e. <math display="inline">n_s=3</math> for 3D problems) and <math display="inline">{}^{n+1}\sigma _{ij}</math> are the Cauchy stresses in <math display="inline">{}^{n+1}V</math>  that are split in the deviatoric (<math display="inline">s_{ij}</math>) and pressure (<math display="inline">p</math>) components as
66
67
<span id="eq-2"></span>
68
{| class="formulaSCP" style="width: 100%; text-align: left;" 
69
|-
70
| 
71
{| style="text-align: left; margin:auto;" 
72
|-
73
| style="text-align: center;" | <math>{}^{n+1}\sigma _{ij} = {}^{n+1}s_{ij}+ {}^{n+1}p \delta _{ij} </math>
74
|}
75
| style="width: 5px;text-align: right;" | (2)
76
|}
77
78
where <math display="inline">\delta _{ij}</math> is the Kronecker delta. The pressure is assumed to be positive for a tension state. Summation of terms with repeated indices  is assumed in Eq.([[#eq-1|1]]) and in  the following.
79
80
'''Remark 1.''' The term <math display="inline">\frac{Dv_i}{Dt}</math> in Eq.([[#eq-1|1]]) is the ''material derivative'' of the velocity. This term is typically computed in a Lagrangian framework as
81
82
<span id="eq-3"></span>
83
{| class="formulaSCP" style="width: 100%; text-align: left;" 
84
|-
85
| 
86
{| style="text-align: left; margin:auto;" 
87
|-
88
| style="text-align: center;" | <math>\frac{Dv_i}{Dt} = \frac{ {}^{n+1}v_i- {}^n v_i }{\Delta t}     </math>
89
|}
90
| style="width: 5px;text-align: right;" | (3)
91
|}
92
93
with
94
95
<span id="eq-4"></span>
96
{| class="formulaSCP" style="width: 100%; text-align: left;" 
97
|-
98
| 
99
{| style="text-align: left; margin:auto;" 
100
|-
101
| style="text-align: center;" | <math>{}^{n+1} v_i:=v_i ({}^{n+1} {x},{}^{n+1} t)\quad ,\quad {}^n v_i:=v_i ({}^n{x},{}^n t)       </math>
102
|}
103
| style="width: 5px;text-align: right;" | (4)
104
|}
105
106
where <math display="inline">{}^n v_i</math> is the value of the velocity  at the material point that has the position <math display="inline">{}^n{x}</math> at time <math display="inline">t={}^nt</math> and <math display="inline">{x}</math> is the coordinates vector in a fixed Cartesian  system <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-48'></span>[[#cite-3|[3,7,48]]].
107
108
===2.2 Constitutive equations===
109
110
===''Newtonian fluids''===
111
112
The relationship between the deviatoric Cauchy stresses <math display="inline">{}^{n+1}s_{ij}</math> and the  rates of deformation <math display="inline">\varepsilon _{ij}</math> has the standard form for a Newtonian fluid,
113
114
<span id="eq-5"></span>
115
{| class="formulaSCP" style="width: 100%; text-align: left;" 
116
|-
117
| 
118
{| style="text-align: left; margin:auto;" 
119
|-
120
| style="text-align: center;" | <math>{}^{n+1}s_{ij} = 2\mu \varepsilon '_{ij} \quad \hbox{with}\quad  \varepsilon '_{ij}= \varepsilon _{ij} - {1 \over 3}\varepsilon _v \delta _{ij}\quad \hbox{and}\quad  \varepsilon _{ij} = {1 \over 2}\left({\partial {}^{n+1}v_i \over \partial {}^{n+1}x_j} +  {\partial {}^{n+1}v_j \over \partial {}^{n+1}x_i}  \right) </math>
121
|}
122
| style="width: 5px;text-align: right;" | (5)
123
|}
124
125
where <math display="inline">\mu </math> is the viscosity, <math display="inline">\varepsilon '_{ij}</math> are the deviatoric rates of deformation and <math display="inline">\varepsilon _v</math> is the volumetric strain rate defined as <math display="inline">\varepsilon _v= \varepsilon _{ii}</math>.
126
127
The above constitutive equation  can be written in matrix form as (for 2D problems)
128
129
<span id="eq-6"></span>
130
{| class="formulaSCP" style="width: 100%; text-align: left;" 
131
|-
132
| 
133
{| style="text-align: left; margin:auto;" 
134
|-
135
| style="text-align: center;" | <math>{}^{n+1}{s} = {D} {\boldsymbol \varepsilon } \quad \hbox{with}\quad   {D}=2 \mu \left[                      \begin{array}{cccccc}2/3 & -1/3 & 0 \\                         -1/3 &2/3  &0  \\                       0 & 0 &1/2\\                      \end{array} \right]   </math>
136
|}
137
| style="width: 5px;text-align: right;" | (6)
138
|}
139
140
In Eq.([[#eq-6|6]]) <math display="inline">{s} = [s_{11},s_{22},s_{12}]^T</math> and <math display="inline">{\boldsymbol \varepsilon } = [\varepsilon _{11},\varepsilon _{22},\varepsilon _{12}]^T</math> are the deviatoric Cauchy stress vector and the strain rate vector, respectively.
141
142
===''Hypo-elastic solids''===
143
144
We will assume a hipo-elastic solid with a constitutive equation of the form
145
146
<span id="eq-7"></span>
147
{| class="formulaSCP" style="width: 100%; text-align: left;" 
148
|-
149
| 
150
{| style="text-align: left; margin:auto;" 
151
|-
152
| style="text-align: center;" | <math>\sigma _{ij}^{\nabla ^J} = 2\bar \mu \varepsilon _{ij} + \lambda \varepsilon _v \delta _{ij}  </math>
153
|}
154
| style="width: 5px;text-align: right;" | (7)
155
|}
156
157
where <math display="inline">\bar \mu </math> and <math display="inline">\lambda </math> are the Lame constants and <math display="inline">\sigma _{ij}^{\nabla ^J}</math> is the Jaumann rate of the  Cauchy stresses <span id='citeF-3'></span>[[#cite-3|[3]]].
158
159
Tensor <math display="inline">{\boldsymbol \sigma }^{\nabla ^J}</math> is approximated in this work as
160
161
<span id="eq-8"></span>
162
{| class="formulaSCP" style="width: 100%; text-align: left;" 
163
|-
164
| 
165
{| style="text-align: left; margin:auto;" 
166
|-
167
| style="text-align: center;" | <math>{\boldsymbol \sigma }^{\nabla ^J}\simeq \frac{1}{J} L ({\boldsymbol \tau })  </math>
168
|}
169
| style="width: 5px;text-align: right;" | (8)
170
|}
171
172
where <math display="inline">J</math> is the determinant of the deformation tensor <math display="inline">{F}</math> <math display="inline">\left(F_{ij} ={\partial {}^{n+1}x_i \over \partial ^{n}x_j}\right)</math>, and <math display="inline">L ({\boldsymbol \tau }) </math> is the Lie derivative of the Kirchhoff stress tensor <math display="inline">\tau _{ij}</math> <span id='citeF-3'></span>[[#cite-3|[3]]].
173
174
Combining Eqs.([[#eq-7|7]]) and ([[#eq-8|8]]) gives the following approximation for <math display="inline">L ({\boldsymbol \tau })</math>
175
176
<span id="eq-9"></span>
177
{| class="formulaSCP" style="width: 100%; text-align: left;" 
178
|-
179
| 
180
{| style="text-align: left; margin:auto;" 
181
|-
182
| style="text-align: center;" | <math>L(\tau _{ij}) \simeq J (2\bar \mu \varepsilon _{ij} + \lambda \varepsilon _v \delta _{ij})  </math>
183
|}
184
| style="width: 5px;text-align: right;" | (9)
185
|}
186
187
'''Remark 2.''' Assumption ([[#eq-8|8]]) is equivalent to accepting that the Truesdell and Jaumann rates of the Cauchy stresses are identical. A more accurate assumption for solving the same problem can be found in <span id='citeF-10'></span>[[#cite-10|[10]]].
188
189
Using standard transformations of continuum mechanics gives <span id='citeF-3'></span>[[#cite-3|[3]]]
190
191
<span id="eq-10"></span>
192
{| class="formulaSCP" style="width: 100%; text-align: left;" 
193
|-
194
| 
195
{| style="text-align: left; margin:auto;" 
196
|-
197
| style="text-align: center;" | <math>{L} ({\boldsymbol \tau }) =  {F}\dot {S} {F}^T = \frac{1}{\Delta t} {F} (\Delta {S}) {F}^T = \frac{J}{\Delta t}\Delta{\boldsymbol \sigma }  </math>
198
|}
199
| style="width: 5px;text-align: right;" | (10)
200
|}
201
202
where <math display="inline">{S}</math> are the 2nd Piola-Kirchhoff stress tensor, the upper dot denotes the time derivative and <math display="inline"> \Delta{\boldsymbol \sigma }</math> is the increment of the Cauchy stress tensor.
203
204
Substituting Eq.([[#eq-10|10]]) into ([[#eq-9|9]])  gives after grouping terms
205
206
<span id="eq-11"></span>
207
{| class="formulaSCP" style="width: 100%; text-align: left;" 
208
|-
209
| 
210
{| style="text-align: left; margin:auto;" 
211
|-
212
| style="text-align: center;" | <math>{}^{n+1}\sigma _{ij}= {}^{n}\sigma _{ij} + \Delta \sigma _{ij}  </math>
213
|}
214
| style="width: 5px;text-align: right;" | (11)
215
|}
216
217
where
218
219
<span id="eq-12"></span>
220
{| class="formulaSCP" style="width: 100%; text-align: left;" 
221
|-
222
| 
223
{| style="text-align: left; margin:auto;" 
224
|-
225
| style="text-align: center;" | <math>\Delta \sigma _{ij} = \Delta s_{ij}+ \Delta p \delta _{ij}  </math>
226
|}
227
| style="width: 5px;text-align: right;" | (12)
228
|}
229
230
with
231
232
<span id="eq-13"></span>
233
{| class="formulaSCP" style="width: 100%; text-align: left;" 
234
|-
235
| 
236
{| style="text-align: left; margin:auto;" 
237
|-
238
| style="text-align: center;" | <math>\Delta s_{ij} = 2\bar \mu \Delta t \varepsilon '_{ij}\quad \hbox{and}\quad \Delta p = \bar \kappa \Delta t \varepsilon _v  </math>
239
|}
240
| style="width: 5px;text-align: right;" | (13)
241
|}
242
243
where
244
245
<span id="eq-14"></span>
246
{| class="formulaSCP" style="width: 100%; text-align: left;" 
247
|-
248
| 
249
{| style="text-align: left; margin:auto;" 
250
|-
251
| style="text-align: center;" | <math>\bar \kappa = \frac{2}{3} \bar \mu + \lambda    </math>
252
|}
253
| style="width: 5px;text-align: right;" | (14)
254
|}
255
256
is the bulk modulus of the solid material.
257
258
The Cauchy stresses <math display="inline">{}^{n}{\boldsymbol \sigma }</math> in Eq.([[#eq-12|12]]) are computed on the updated configuration. These stresses can be obtained from the 2nd Piola-Kirchhoff stresses <math display="inline">{}^n{S}</math> on the current configuration  as
259
260
<span id="eq-15"></span>
261
{| class="formulaSCP" style="width: 100%; text-align: left;" 
262
|-
263
| 
264
{| style="text-align: left; margin:auto;" 
265
|-
266
| style="text-align: center;" | <math>{}^{n}{\boldsymbol \sigma } = \frac{1}{J} {F}  ({}^n{S}){F}^T   </math>
267
|}
268
| style="width: 5px;text-align: right;" | (15)
269
|}
270
271
'''Remark 3.''' The relationship between the deviatoric Cauchy stress increments and the rates of deformation for a solid material can be written in the same matrix form as Eq.(6) substituting in matrix <math display="inline">D</math> the viscosity <math display="inline">\mu </math> by <math display="inline">\bar \mu \Delta t </math>.
272
273
===2.3 Mass conservation equation===
274
275
The  mass conservation equation can be written for both fluid and solids as <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-48'></span>[[#cite-3|[3,7,48]]]
276
277
<span id="eq-16"></span>
278
{| class="formulaSCP" style="width: 100%; text-align: left;" 
279
|-
280
| 
281
{| style="text-align: left; margin:auto;" 
282
|-
283
| style="text-align: center;" | <math>r_v =0  \quad  \hbox{in } {}^{n+1}V  </math>
284
|}
285
| style="width: 5px;text-align: right;" | (16)
286
|}
287
288
with
289
290
<span id="eq-17"></span>
291
{| class="formulaSCP" style="width: 100%; text-align: left;" 
292
|-
293
| 
294
{| style="text-align: left; margin:auto;" 
295
|-
296
| style="text-align: center;" | <math>r_v:=-\frac{1}{\kappa }\frac{Dp}{Dt}+\varepsilon _v </math>
297
|}
298
| style="width: 5px;text-align: right;" | (17)
299
|}
300
301
For a quasi-incompressible fluid material <math display="inline">\kappa =\rho c^2</math> where <math display="inline">c</math> is the speed of sound in the fluid. For a solid material, <math display="inline">\kappa = \bar \kappa </math> with <math display="inline">\bar \kappa </math> given by Eq.([[#eq-14|14]]). Eqs.([[#2.3 Mass conservation equation|2.3]]) define the constitutive equation for the pressure for both solids and quasi-incompressible fluids as <math display="inline">\Delta p = \kappa  \Delta t \varepsilon _v</math>. This is consisted with the expression obtained in Eq.([[#eq-13|13]]) for a hypo-elastic solid material.
302
303
For a fully incompressible material <math display="inline">\kappa =\infty </math> and Eq.([[#eq-16|16]]) simplifies to the standard form, <math display="inline">\varepsilon _v =0</math>. In this case the pressure increment can not be obtained from Eqs.([[#2.3 Mass conservation equation|2.3]]) and then the pressure becomes an independent variable <span id='citeF-7'></span><span id='citeF-48'></span>[[#cite-7|[7,48]]]. In our work we will retain the quasi-incompressible form of <math display="inline">r_v</math>  for convenience.
304
305
'''Remark 4.''' Note that for Newtonian fluids the deviatoric stresses are directly related to the rates of deformation by Eq.([[#eq-5|5]]), whereas for solids the ''deviatoric stress increments'' are related to the rates of deformation by Eq.([[#eq-13|13]]). On the other hand, the relationship between the pressure increment and the volumetric strain rate has the form of Eq.([[#eq-13|13]]) for both fluids and solids, with the adequate definition of the bulk modulus for each case.
306
307
===2.4 Thermal balance===
308
309
The thermal balance equation is written in a Lagrangian framework as
310
311
<span id="eq-18"></span>
312
{| class="formulaSCP" style="width: 100%; text-align: left;" 
313
|-
314
| 
315
{| style="text-align: left; margin:auto;" 
316
|-
317
| style="text-align: center;" | <math>\rho c \frac{DT}{Dt} - {\partial  \over \partial {}^{n+1} x_i} \left(k {\partial T \over \partial {}^{n+1} x_i}\right)+ {}^{n+1} Q =0  \quad ,\quad  i,j=1,\cdots , n_s  \quad \hbox{in } {}^{n+1}V  </math>
318
|}
319
| style="width: 5px;text-align: right;" | (18)
320
|}
321
322
where <math display="inline">T</math> is the temperature, <math display="inline">c</math> is the thermal capacity, <math display="inline">k</math> is the heat conductivity and <math display="inline">Q</math> is the external heat source.
323
324
===2.5 Boundary conditions===
325
326
===''Mechanical problem''===
327
328
The boundary conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann (<math display="inline">\Gamma _t</math>) boundaries with <math display="inline">\Gamma  = \Gamma _v \cup \Gamma _t</math> are
329
330
<span id="eq-19"></span>
331
{| class="formulaSCP" style="width: 100%; text-align: left;" 
332
|-
333
| 
334
{| style="text-align: left; margin:auto;" 
335
|-
336
| style="text-align: center;" | <math>{}^{n+1} v_i -{}^{n+1} v_i^p =0 \qquad \hbox{on }{}^{n+1}\Gamma _v </math>
337
|}
338
| style="width: 5px;text-align: right;" | (19)
339
|}
340
341
<span id="eq-20"></span>
342
{| class="formulaSCP" style="width: 100%; text-align: left;" 
343
|-
344
| 
345
{| style="text-align: left; margin:auto;" 
346
|-
347
| style="text-align: center;" | <math>{}^{n+1}\sigma _{ij}n_j - {}^{n+1} t_i^p =0 \quad \hbox{on }{}^{n+1}\Gamma _t \quad , \quad i,j=1,\cdots , n_s  </math>
348
|}
349
| style="width: 5px;text-align: right;" | (20)
350
|}
351
352
where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math> are the prescribed velocities and prescribed tractions on the respective boundaries and <math display="inline">n_j</math> are the components of the unit normal vector to the boundary <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-48'></span>[[#cite-3|[3,7,48]]].
353
354
===''Thermal problem''===
355
356
<span id="eq-21"></span>
357
<span id="eq-22"></span>
358
{| class="formulaSCP" style="width: 100%; text-align: left;" 
359
|-
360
| 
361
{| style="text-align: left; margin:auto;" 
362
|-
363
| style="text-align: center;" | <math>{}^{n+1}T - {}^{n+1}T^p =0 \quad \hbox{on }{}^{n+1}\Gamma _T</math>
364
| style="width: 5px;text-align: right;" | (21)
365
|-
366
| style="text-align: center;" | <math> k {\partial T \over \partial n} + {}^{n+1} q_n^p =0 \quad \hbox{on }{}^{n+1}\Gamma _q  </math>
367
| style="width: 5px;text-align: right;" | (22)
368
|}
369
|}
370
371
where <math display="inline">T^p</math> and <math display="inline">q_n^p</math> are the prescribed temperature and the prescribed normal heat flux at the boundaries <math display="inline">\Gamma _T</math> and <math display="inline">\Gamma _q</math>, respectively and <math display="inline">n</math> is the direction normal to be boundary.
372
373
==3 STABILIZED MASS CONSERVATION EQUATION==
374
375
In this work we will use the following stabilized form of the mass conservation equation obtained using the Finite Calculus (FIC) approach <span id='citeF-38'></span>[[#cite-38|[38]]]
376
377
<span id="eq-23"></span>
378
{| class="formulaSCP" style="width: 100%; text-align: left;" 
379
|-
380
| 
381
{| style="text-align: left; margin:auto;" 
382
|-
383
| style="text-align: center;" | <math>-\frac{1}{\kappa }\frac{Dp}{Dt}+  \varepsilon _v  +\tau {\partial {}^{n+1} \hat r_{m_i} \over \partial {}^{n+1}x_i} =0 \qquad \hbox{in }{}^{n+1}V  </math>
384
|}
385
| style="width: 5px;text-align: right;" | (23)
386
|}
387
388
The last term in Eq.([[#eq-23|23]]) is the so-called stabilization term introduced by the FIC technique; in this equation <math display="inline">\hat r_{m_i}</math> is a ''static momentum term'' defined as
389
390
<span id="eq-24"></span>
391
{| class="formulaSCP" style="width: 100%; text-align: left;" 
392
|-
393
| 
394
{| style="text-align: left; margin:auto;" 
395
|-
396
| style="text-align: center;" | <math>\hat r_{m_i}= {\partial s_{ij} \over \partial x_j} + \frac{\partial p}{\partial x_i} + b_i   </math>
397
|}
398
| style="width: 5px;text-align: right;" | (24)
399
|}
400
401
and <math display="inline">\tau </math> is a ''stabilization parameter''. Eq.([[#eq-23|23]]) is a simplification of a more general FIC-based stabilized expression for quasi-incompressible fluids detailed  in <span id='citeF-38'></span>[[#cite-38|[38]]].
402
403
The stabilization parameter <math display="inline">\tau </math> is typically computed for each element  as
404
405
<span id="eq-25"></span>
406
{| class="formulaSCP" style="width: 100%; text-align: left;" 
407
|-
408
| 
409
{| style="text-align: left; margin:auto;" 
410
|-
411
| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{(l^e)^2} + \frac{2\rho }{\Delta t}  \right)^{-1} </math>
412
|}
413
| style="width: 5px;text-align: right;" | (25)
414
|}
415
416
where <math display="inline">\Delta t</math> is the time step used for the transient solution and <math display="inline">l^e</math> is a characteristic element length computed as <math display="inline">l^e = 2(V^e)^{1/n_s}</math> where <math display="inline">V^e</math> is the element area (for 3-noded triangles) or volume (for 4-noded tetrahedra). For heterogeneous material the values of <math display="inline">\mu </math> and <math display="inline">\rho </math>  are computed at the element center <span id='citeF-38'></span>[[#cite-38|[38]]]. For  a solid material, the viscosity <math display="inline">\mu </math> is substituted by <math display="inline">\bar \mu \Delta t</math> with <math display="inline">\bar \mu </math> defined in Eq.(7).
417
418
We note that <math display="inline">\tau =0</math> for “compressible” type problems not requiring stabilization. In these cases Eq.([[#eq-23|23]]) recovers the standard form for the conservation of mass of the infinitesimal theory (see Eqs.([[#2.3 Mass conservation equation|2.3]])).
419
420
==4 VARIATIONAL EQUATIONS==
421
422
===4.1 Momentum equations===
423
424
Multiplying Eq.(1) by arbitrary test functions <math display="inline">w_i</math> with dimensions of velocity and integrating over the analysis domain <math display="inline">{}^{n+1}V</math> gives the weighted residual form of the momentum equations as <span id='citeF-46'></span><span id='citeF-48'></span>[[#cite-46|[46,48]]]
425
426
<span id="eq-26"></span>
427
{| class="formulaSCP" style="width: 100%; text-align: left;" 
428
|-
429
| 
430
{| style="text-align: left; margin:auto;" 
431
|-
432
| style="text-align: center;" | <math>\int _{{}^{n+1}V} w_i \left(\rho \frac{Dv_i}{Dt}- {\partial {}^{n+1}\sigma _{ij} \over \partial {}^{n+1} x_j}- {}^{n+1} b_i\right)dV =0 </math>
433
|}
434
| style="width: 5px;text-align: right;" | (26)
435
|}
436
437
Integrating by parts  the term involving <math display="inline">\sigma _{ij}</math> and using the Neumann boundary conditions ([[#eq-20|20]]) yields the weak variational form of the momentum equations as
438
439
<span id="eq-27"></span>
440
{| class="formulaSCP" style="width: 100%; text-align: left;" 
441
|-
442
| 
443
{| style="text-align: left; margin:auto;" 
444
|-
445
| style="text-align: center;" | <math>\int _{{}^{n+1}V} w_i \rho \frac{Dv_i}{Dt} dV + \int _{{}^{n+1}V} \delta \varepsilon _{ij} {}^{n+1} \sigma _{ij} dV - \int _{{}^{n+1}V} w_i {}^{n+1} b_i dV - \int _{{}^{n+1}\Gamma _t} w_i {}^{n+1} t_i^p d\Gamma =0 </math>
446
|}
447
| style="width: 5px;text-align: right;" | (27)
448
|}
449
450
where <math display="inline">\delta \varepsilon _{ij}= {\partial w_i \over \partial {}^{n+1}x_j}+{\partial w_j \over \partial {}^{n+1}x_i}</math> is an arbitrary (virtual) strain rate field. Eq.([[#eq-27|27]]) is the standard form of the ''principle of virtual power''  <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-48'></span>[[#cite-3|[3,7,48]]].
451
452
Substituting  the stresses from Eq.([[#eq-2|2]]) into ([[#eq-27|27]]) gives the following expression (in matrix form)
453
454
<span id="eq-28"></span>
455
{| class="formulaSCP" style="width: 100%; text-align: left;" 
456
|-
457
| 
458
{| style="text-align: left; margin:auto;" 
459
|-
460
| style="text-align: center;" | <math>\int _{{}^{n+1}V} {w}^T \rho \frac{D\hbox{v}}{Dt} dV + \!\int _{{}^{n+1}V} \delta {\boldsymbol \varepsilon }^T {}^{n+1}{s} dV  + \!\int _{{}^{n+1}V} \delta {\boldsymbol \varepsilon }^T {m} {}^{n+1}{p} dV -</math>
461
|-
462
| style="text-align: center;" | <math>  -\int _{{}^{n+1}V} {w}^T {}^{n+1}{b}dV - \!\int _{{}^{n+1}\Gamma _t} \hbox{w}^T {}^{n+1}\hbox{t}^p d\Gamma =0  </math>
463
|}
464
| style="width: 5px;text-align: right;" | (28)
465
|}
466
467
In Eq.([[#eq-28|28]]) <math display="inline">{w},{v}</math> and <math display="inline">\delta {\boldsymbol \varepsilon }</math> are vectors containing the test functions, the velocities and the virtual strain rates respectively; <math display="inline">{b}</math> and <math display="inline">\hbox{t}^p</math> are the body force and surface traction vectors, respectively and <math display="inline">{m}</math> is an auxiliary vector. These vectors are defined as (for 2D problems)
468
469
<span id="eq-29"></span>
470
{| class="formulaSCP" style="width: 100%; text-align: left;" 
471
|-
472
| 
473
{| style="text-align: left; margin:auto;" 
474
|-
475
| style="text-align: center;" | <math>\begin{array}{l}\hbox{w} =[w_1,w_2]^T\quad ,\quad \hbox{v} =[v_1,v_2]^T \quad ,\quad \hbox{b} =[b_1,b_2]^T \quad ,\quad  \hbox{t}^p =[t_1^p,t_2^p]^T\\[.5cm]  \delta {\boldsymbol \varepsilon }= [\delta \varepsilon _{11},\delta \varepsilon _{22}, \delta \varepsilon _{12}]^T \quad ,\quad  {m}=[1,1,0]^T  \end{array} </math>
476
|}
477
| style="width: 5px;text-align: right;" | (29)
478
|}
479
480
===4.2 Mass conservation equation===
481
482
We multiply Eq.([[#eq-23|23]]) by arbitrary (continuous) test functions <math display="inline">q</math> (with dimensions of pressure) defined over the analysis domain. Integrating over <math display="inline">{}^{n+1}V</math> gives
483
484
<span id="eq-30"></span>
485
{| class="formulaSCP" style="width: 100%; text-align: left;" 
486
|-
487
| 
488
{| style="text-align: left; margin:auto;" 
489
|-
490
| style="text-align: center;" | <math>\int _{{}^{n+1}V} -\frac{q}{\kappa } \frac{Dp}{Dt}dV +  \int _{{}^{n+1}V} q \varepsilon _v dV + \int _{{}^{n+1}V} q \tau {\partial {}^{n+1}\hat r_{m_i} \over \partial {}^{n+1}x_i} dV =0  </math>
491
|}
492
| style="width: 5px;text-align: right;" | (30)
493
|}
494
495
The third integral in Eq.([[#eq-30|30]]) corresponding to the stabilization term is zero in the compressible parts of the continuum.
496
497
Integrating by parts the last integral in Eq.([[#eq-30|30]]) and using ([[#eq-24|24]]) gives after some transformations <span id='citeF-38'></span>[[#cite-38|[38]]]
498
499
<span id="eq-31"></span>
500
{| class="formulaSCP" style="width: 100%; text-align: left;" 
501
|-
502
| 
503
{| style="text-align: left; margin:auto;" 
504
|-
505
| style="text-align: center;" | <math>\begin{array}{r}\displaystyle \int _{{}^{n+1}V} \frac{q}{\kappa } \frac{Dp}{Dt}dV - \int _{{}^{n+1}V}q \varepsilon _v dV + \int _{{}^{n+1}V} \tau {\partial q \over \partial {}^{n+1}x_i} \left({\partial {}^{n+1}s_{ij} \over \partial {}^{n+1}x_j} + {\partial {}^{n+1}p \over \partial {}^{n+1}x_i} + {}^{n+1}b_i\right)dV\\ \displaystyle - \int _{{}^{n+1}\Gamma _t} q \tau \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} \left({}^{n+1}s_n +{}^{n+1} p - {}^{n+1}t_n\right)\right]d\Gamma =0 \end{array}   </math>
506
|}
507
| style="width: 5px;text-align: right;" | (31)
508
|}
509
510
where <math display="inline">v_n</math> and <math display="inline">t_n</math> are the velocity and the traction force normal to the boundary <math display="inline">\Gamma _t</math>, respectively and <math display="inline">s_n = s_{ij} n_j</math>. The characteristic length <math display="inline">h_n</math> in Eq.([[#eq-31|31]]) is taken as <math display="inline">h_n = \frac{l^e}{2}</math> in practice.
511
512
===4.3 Thermal balance equation===
513
514
Application of the standard weighted residual method to the thermal balance equations ([[#eq-18|18]]), ([[#eq-21|21]]) and ([[#eq-22|22]]) leads, after standard operations, to <span id='citeF-46'></span><span id='citeF-48'></span>[[#cite-46|[46,48]]]
515
516
<span id="eq-32"></span>
517
{| class="formulaSCP" style="width: 100%; text-align: left;" 
518
|-
519
| 
520
{| style="text-align: left; margin:auto;" 
521
|-
522
| style="text-align: center;" | <math>\int _{{}^{n+1}V} \hat w \rho c \frac{DT}{Dt}dV + \int _{{}^{n+1}V} {\partial \hat w  \over \partial {}^{n+1}x_i}  k {\partial T \over \partial {}^{n+1}x_i}dV - \int _{{}^{n+1}V} \hat w {}^{n+1} Q dV+ \int _{{}^{n+1}\Gamma _q} \hat w {}^{n+1} q_n^p d\Gamma =0 </math>
523
|}
524
| style="width: 5px;text-align: right;" | (32)
525
|}
526
527
where <math display="inline">\hat w</math> are the space weighting functions for the temperature.
528
529
==5 FEM DISCRETIZATION==
530
531
We discretize  the analysis domain into finite elements with <math display="inline">n</math> nodes in the standard manner leading to a mesh with a total number of <math display="inline">N_e</math> elements and <math display="inline">N</math> nodes <span id='citeF-33'></span><span id='citeF-46'></span>[[#cite-33|[33,46]]]. In our work we will choose simple 3-noded linear triangles (<math display="inline">n=3</math>) for 2D problems and 4-noded tetrahedra (<math display="inline">n=4</math>) for 3D problems with local linear shape functions <math display="inline">N_i^e</math> defined for each node <math display="inline">i</math> (<math display="inline">i=1,\cdots , n</math>) of element <math display="inline">e</math> <span id='citeF-33'></span><span id='citeF-47'></span>[[#cite-33|[33,47]]]. The velocity components, the  pressure and the temperature are interpolated over the mesh in terms of their nodal values in the same manner using the  ''global linear shape functions'' <math display="inline">N_j</math> spanning over the elements sharing node <math display="inline">j</math> (<math display="inline">j=1,\cdots , N</math>) <span id='citeF-33'></span><span id='citeF-46'></span>[[#cite-33|[33,46]]]. In matrix form we have (for 2D problems)
532
533
<span id="eq-33"></span>
534
{| class="formulaSCP" style="width: 100%; text-align: left;" 
535
|-
536
| 
537
{| style="text-align: left; margin:auto;" 
538
|-
539
| style="text-align: center;" | <math>{v} = {N}_v\bar {v}  \quad ,\quad p = {N}_p\bar {p}\quad  , \quad T = {N}_T \bar {T}  </math>
540
|}
541
| style="width: 5px;text-align: right;" | (33)
542
|}
543
544
where
545
546
<span id="eq-34"></span>
547
{| class="formulaSCP" style="width: 100%; text-align: left;" 
548
|-
549
| 
550
{| style="text-align: left; margin:auto;" 
551
|-
552
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \bar {v} = \left\{\begin{matrix}\bar{v}^1\\\bar{v}^2\\\vdots \\ \bar{v}^N\end{matrix}  \right\}\quad \hbox{with}~ \bar{v}^i = \left\{\begin{matrix}\bar{v}^i_1\\\bar{v}^i_2\end{matrix}  \right\}~, ~ \bar {p} =\left\{\begin{matrix}\bar{p}^1\\\bar{p}^2\\\vdots \\ \bar{p}^{N}\end{matrix}  \right\}~, ~ \bar {T}= \left\{\begin{matrix}\bar{T}_1\\\bar{T}_2\\\vdots \\ \bar{T}^N\end{matrix}  \right\}\\[.5cm] \displaystyle {N}_v = [{N}_1, {N}_2,\cdots , {N}_N ]^T ~,~ {N}_p = {N}_T = [{N}_1, { N}_2,\cdots , {N}_N ] \end{array} </math>
553
|}
554
| style="width: 5px;text-align: right;" | (34)
555
|}
556
557
with <math display="inline">{N}_j = N_j {I}_2</math> where <math display="inline">{I}_2</math> is the <math display="inline">2\times 2</math> unit matrix.
558
559
In Eq.([[#eq-34|34]]) vectors <math display="inline">\bar{v}</math>, <math display="inline">\bar {p}</math> and <math display="inline">\bar {T}</math>  contain the nodal velocities, the nodal pressures and the nodal temperatures for the whole mesh, respectively and the upperindex denotes the nodal value for each vector or scalar magnitude.
560
561
Substituting Eqs.([[#eq-33|33]]) into Eqs.([[#eq-27|27]]), ([[#eq-31|31]]) and ([[#eq-32|32]]) and choosing a Galerkin formulation with <math display="inline">w_i = q =  \hat w_i =N_i</math> leads to the following system of algebraic equations
562
563
<span id="eq-35"></span>
564
{| class="formulaSCP" style="width: 100%; text-align: left;" 
565
|-
566
| 
567
{| style="text-align: left; margin:auto;" 
568
|-
569
| style="text-align: center;" | <math>{M}_0 {\dot{\bar{v}}} + {g} +{Q}{}^{n+1}\bar {p}- {f}_v={0} </math>
570
|}
571
| style="width: 5px;text-align: right;" | (35)
572
|}
573
574
<span id="eq-36"></span>
575
{| class="formulaSCP" style="width: 100%; text-align: left;" 
576
|-
577
| 
578
{| style="text-align: left; margin:auto;" 
579
|-
580
| style="text-align: center;" | <math>\hbox{M}_1 {\dot{\bar{p}}}-{Q}^T {}^{n+1}\bar{v} + ({L}+{M}_b){}^{n+1} \bar {p}- {f}_p={0} </math>
581
|}
582
| style="width: 5px;text-align: right;" | (36)
583
|}
584
585
<span id="eq-37"></span>
586
{| class="formulaSCP" style="width: 100%; text-align: left;" 
587
|-
588
| 
589
{| style="text-align: left; margin:auto;" 
590
|-
591
| style="text-align: center;" | <math>{C} {\dot{\bar{T}}}+\hat {L} {}^{n+1}\bar{T} -  {f}_T={0}  </math>
592
|}
593
| style="width: 5px;text-align: right;" | (37)
594
|}
595
596
where <math display="inline">{\dot{\bar{a}}}</math>  denotes the material time derivative of the components of  a vector <math display="inline">{a}</math>. The element form of the matrices and vectors in Eqs.([[#5 FEM DISCRETIZATION|5]]) is given in Box 1.
597
598
We note that Eqs.([[#5 FEM DISCRETIZATION|5]]) involve the geometry at the updated configuration (<math display="inline">{}^{n+1}V</math>). This geometry is unknown; hence the solution of Eqs.([[#5 FEM DISCRETIZATION|5]]) has to be found iteratively. The iterative solution scheme chosen in this work is presented in the next section.
599
600
<br/>
601
602
'''Remark 5.''' The boundary terms of vector <math display="inline">{f}_p</math> can be incorporated into the matrices of Eq.([[#eq-36|36]]). This leads to a non symmetrical set of equations. These boundary terms are computed here iteratively within the incremental solution scheme.
603
604
'''Remark 6.''' The presence of matrix <math display="inline">{M}_b</math> in Eq.([[#eq-36|36]]) allows us to compute the pressure without the need of prescribing its value at the free surface <span id='citeF-38'></span>[[#cite-38|[38]]].
605
606
==6 TRANSIENT SOLUTION OF THE DISCRETIZED EQUATIONS==
607
608
Eqs.([[#5 FEM DISCRETIZATION|5]]) are solved in time with an implicit Newton-Raphson type iterative scheme <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-47'></span>[[#cite-3|[3,7,47]]]. The basic steps within a time increment <math display="inline">[n,n+1]</math> are:<br/>
609
610
* Initialize variables: <math display="inline">({}^{n+1}\hbox{x}^1,{}^{n+1}\bar{v}^1,{}^{n+1}{s}^1,{}^{n+1}\bar{p}^1,  {}^{n+1}\bar{T}^1)\equiv  \left\{{}^{n}{x},{}^{n}\bar{v},{}^{n}{s},{}^{n}\bar{p},{}^{n}\bar{T}\right\}</math>.<p> In the following <math display="inline">(\cdot )^i</math> denotes a value computed at the <math display="inline">i</math>th iteration. </p>
611
* Iteration loop: <math display="inline">i=1,\cdots , N_{iter}</math>.<p> For each iteration. 
612
613
</p>
614
615
===''Step 1. Compute the initial Cauchy stresses <math>{}^{n}{\boldsymbol \sigma }^i</math> for solid elements''===
616
617
<span id="eq-38"></span>
618
{| class="formulaSCP" style="width: 100%; text-align: left;" 
619
|-
620
| 
621
{| style="text-align: left; margin:auto;" 
622
|-
623
| style="text-align: center;" | <math>{}^{n} {\boldsymbol \sigma }^i =  \frac{1}{J^i}{F}^i [{}^{n}{S}] {F}^{iT} \quad \hbox{with } F_{kl}^i = {\partial {}^{n+1}x_k^i \over \partial {}^{n} x_l} \quad \hbox{and } {J^i} = \vert {F}^i\vert   </math>
624
|}
625
| style="width: 5px;text-align: right;" | (38)
626
|}
627
628
with
629
630
<span id="eq-39"></span>
631
{| class="formulaSCP" style="width: 100%; text-align: left;" 
632
|-
633
| 
634
{| style="text-align: left; margin:auto;" 
635
|-
636
| style="text-align: center;" | <math>{}^{n}{S} = {}^{n} {\boldsymbol \sigma }\quad \hbox{and}\quad {}^{n}{\sigma }_{ij} = {}^{n}s_{ij} + {}^{n}p \delta _{ij}  </math>
637
|}
638
| style="width: 5px;text-align: right;" | (39)
639
|}
640
641
===''Step 2. Compute the nodal velocity increments <math>\Delta \bar {v}</math>''===
642
643
From Eq.([[#5 FEM DISCRETIZATION|5]]a)
644
645
<span id="eq-40"></span>
646
{| class="formulaSCP" style="width: 100%; text-align: left;" 
647
|-
648
| 
649
{| style="text-align: left; margin:auto;" 
650
|-
651
| style="text-align: center;" | <math>{H}_v^i \Delta \bar {v} = - {}^{n+1}\bar{\hbox{r}}_m^i \rightarrow  \Delta \bar{v} </math>
652
|}
653
| style="width: 5px;text-align: right;" | (40)
654
|}
655
656
with the momentum residual <math display="inline">\bar{r}_m</math> and the iteration matrix <math display="inline">{H}_v</math>  given by
657
658
<span id="eq-41"></span>
659
{| class="formulaSCP" style="width: 100%; text-align: left;" 
660
|-
661
| 
662
{| style="text-align: left; margin:auto;" 
663
|-
664
| style="text-align: center;" | <math>\bar{\hbox{r}}_m = \hbox{M}_0 {\dot{\bar{v}}} + \hbox{g}+\hbox{Q}\bar {p}-\hbox{f}_v\quad , \quad \hbox{H}_v = \frac{1}{\Delta t} \hbox{M}_0 + \hbox{K} + \hbox{K}_v  </math>
665
|}
666
| style="width: 5px;text-align: right;" | (41)
667
|}
668
669
where <math display="inline"> \hbox{K}_{v}</math> is the bulk tangent matrix (see Remark 9) and
670
671
<span id="eq-42"></span>
672
{| class="formulaSCP" style="width: 100%; text-align: left;" 
673
|-
674
| 
675
{| style="text-align: left; margin:auto;" 
676
|-
677
| style="text-align: center;" | <math>\hbox{K}_{ij}^e = \int _{{}^nV^e} \left(\hbox{B}_i^{eT} {\boldsymbol{\cal C}}_T \hbox{B}_j^{e} dV + {G}^T_i [{S} ]  {G}_j \Delta t \right)dV   </math>
678
|}
679
| style="width: 5px;text-align: right;" | (42)
680
|}
681
682
where <math display="inline">{\boldsymbol{\cal C}}_T</math> is the tangent constitutive matrix emanating from the linearization of Eq.([[#eq-35|35]]) (expressed in the current configuration) with respect to the nodal velocities <span id='citeF-10'></span><span id='citeF-37'></span>[[#cite-10|[10,37]]]. The second term in matrix <math display="inline">\hbox{K}^e</math> is explained in  Remark 7.
683
684
===''Step 3. Update the nodal velocities''===
685
686
<span id="eq-43"></span>
687
{| class="formulaSCP" style="width: 100%; text-align: left;" 
688
|-
689
| 
690
{| style="text-align: left; margin:auto;" 
691
|-
692
| style="text-align: center;" | <math>{}^{n+1}\bar{v}^{i+1}= \bar{v}^i + \Delta \bar{v} </math>
693
|}
694
| style="width: 5px;text-align: right;" | (43)
695
|}
696
697
===''Step 4. Compute the nodal pressures''===
698
699
From Eq.([[#eq-36|36]]) we obtain
700
701
<span id="eq-44"></span>
702
{| class="formulaSCP" style="width: 100%; text-align: left;" 
703
|-
704
| 
705
{| style="text-align: left; margin:auto;" 
706
|-
707
| style="text-align: center;" | <math>\hbox{H}_p^i  {}^{n+1}\bar{p}^{i+1}= \frac{1}{\Delta t} \hbox{M}_1\bar{p}^i  + {Q}^T \bar{v}^{i+1}+ \bar{f}^{i}_p  \rightarrow {}^{n+1}\bar{p}^{i+1}   </math>
708
|}
709
| style="width: 5px;text-align: right;" | (44)
710
|}
711
712
with
713
714
<span id="eq-45"></span>
715
{| class="formulaSCP" style="width: 100%; text-align: left;" 
716
|-
717
| 
718
{| style="text-align: left; margin:auto;" 
719
|-
720
| style="text-align: center;" | <math>\hbox{H}_p = \frac{1}{\Delta t} \hbox{M}_1 + \hbox{L} + \hbox{M}_b </math>
721
|}
722
| style="width: 5px;text-align: right;" | (45)
723
|}
724
725
===''Step 5. Update the nodal coordinates and the deformation gradient''===
726
727
<span id="eq-46"></span>
728
{| class="formulaSCP" style="width: 100%; text-align: left;" 
729
|-
730
| 
731
{| style="text-align: left; margin:auto;" 
732
|-
733
| style="text-align: center;" | <math>{}^{n+1}{x}^{i+1}= {x}^i + \frac{1}{2} (\bar{v}^{i+1} + {}^{n}\bar{v})\Delta t </math>
734
|}
735
| style="width: 5px;text-align: right;" | (46)
736
|}
737
738
{| class="formulaSCP" style="width: 100%; text-align: left;" 
739
|-
740
| 
741
{| style="text-align: left; margin:auto;" 
742
|-
743
| style="text-align: center;" | <math>{}^{n+1}F^{i+1}_{ij} = {\partial {}^{n+1}{x}^{i+1}_i \over \partial {}^{n}{x_j}}</math>
744
|}
745
|}
746
747
A more accurate expression for computing <math display="inline">{}^{n+1}{x}^{i+1}</math> can be used involving the  nodal accelerations <span id='citeF-37'></span>[[#cite-37|[37]]].
748
749
===''Step 6.Compute the deviatoric Cauchy stresses''===
750
751
''Fluids:''
752
753
<span id="eq-47"></span>
754
{| class="formulaSCP" style="width: 100%; text-align: left;" 
755
|-
756
| 
757
{| style="text-align: left; margin:auto;" 
758
|-
759
| style="text-align: center;" | <math>{}^{n+1}{s}^{i+1} = {D}_T\hbox{B}\bar {v}^{i+1}  </math>
760
|}
761
| style="width: 5px;text-align: right;" | (47)
762
|}
763
764
''Solids:''
765
766
<span id="eq-48"></span>
767
{| class="formulaSCP" style="width: 100%; text-align: left;" 
768
|-
769
| 
770
{| style="text-align: left; margin:auto;" 
771
|-
772
| style="text-align: center;" | <math>\Delta  {s}= {D}_T  \hbox{B} \Delta \bar {v}  </math>
773
|}
774
| style="width: 5px;text-align: right;" | (48)
775
|}
776
777
where <math display="inline">{D}_T</math> is the tangent constitutive matrix for the deviatoric Cauchy stresses. For elastic solids and Newtonian fluids <math display="inline">{D}_T = {D}</math> where <math display="inline">{D}</math> is given by Eq.([[#eq-6|6]]) and <math display="inline">\mu =\bar \mu \Delta t</math> for solids.
778
779
For non-linear continua, the adequate expression for <math display="inline">{D}_T</math> has to be used <span id='citeF-3'></span><span id='citeF-10'></span>[[#cite-3|[3,10]]].
780
781
===''Step 7.Compute the  stresses''===
782
783
''Solids:''<br/>
784
785
Cauchy stresses:
786
787
<span id="eq-49"></span>
788
{| class="formulaSCP" style="width: 100%; text-align: left;" 
789
|-
790
| 
791
{| style="text-align: left; margin:auto;" 
792
|-
793
| style="text-align: center;" | <math>{}^{n+1}\sigma _{ij}^{i+1} = {}^{n}\sigma _{ij}^{i} + \Delta  {s}_{ij} + \Delta \bar{p}\delta _{ij}  </math>
794
|}
795
| style="width: 5px;text-align: right;" | (49)
796
|}
797
798
with
799
800
<span id="eq-50"></span>
801
{| class="formulaSCP" style="width: 100%; text-align: left;" 
802
|-
803
| 
804
{| style="text-align: left; margin:auto;" 
805
|-
806
| style="text-align: center;" | <math>\Delta \bar{p}= {}^{n+1}\bar{p}^{i+1} - {}^{n}\bar {p}  </math>
807
|}
808
| style="width: 5px;text-align: right;" | (50)
809
|}
810
811
2d Piola Kirchhoff stresses:
812
813
<span id="eq-51"></span>
814
{| class="formulaSCP" style="width: 100%; text-align: left;" 
815
|-
816
| 
817
{| style="text-align: left; margin:auto;" 
818
|-
819
| style="text-align: center;" | <math>{}^{n+1}{S}^{i+1} = J {F}^{-1}[{}^{n+1}{\boldsymbol \sigma }^{i+1}] {F}^{-T} ~~\hbox{with } ~~{F} \equiv {}^{n+1}{F}^{i+1}  </math>
820
|}
821
| style="width: 5px;text-align: right;" | (51)
822
|}
823
824
''Fluids:''
825
826
<span id="eq-52"></span>
827
{| class="formulaSCP" style="width: 100%; text-align: left;" 
828
|-
829
| 
830
{| style="text-align: left; margin:auto;" 
831
|-
832
| style="text-align: center;" | <math>{}^{n+1}\sigma _{ij}^{i+1} = {}^{n+1}{s}^{i+1}_{ij} + {}^{n+1}{p}^{i+1}\delta _{ij}  </math>
833
|}
834
| style="width: 5px;text-align: right;" | (52)
835
|}
836
837
===''Step 8.Compute the nodal temperatures''===
838
839
<span id="eq-53"></span>
840
{| class="formulaSCP" style="width: 100%; text-align: left;" 
841
|-
842
| 
843
{| style="text-align: left; margin:auto;" 
844
|-
845
| style="text-align: center;" | <math>\left[\frac{1}{\Delta t}{C} + \hat{L}   \right]\Delta \bar{T}=- {}^{n+1} \bar{r}^i_T \quad ,\quad  {}^{n+1}\bar{T}^{i+1}= {}^{n}\bar{T}^i + \Delta \bar{T} </math>
846
|}
847
| style="width: 5px;text-align: right;" | (53)
848
|}
849
850
with
851
852
<span id="eq-54"></span>
853
{| class="formulaSCP" style="width: 100%; text-align: left;" 
854
|-
855
| 
856
{| style="text-align: left; margin:auto;" 
857
|-
858
| style="text-align: center;" | <math>\bar {r}_T= {C} {\dot{\bar{T}}} + \hat{L} \bar{T} - {f}_T </math>
859
|}
860
| style="width: 5px;text-align: right;" | (54)
861
|}
862
863
===''Step 9. Check convergence''===
864
865
Verify the following conditions:
866
867
<span id="eq-55"></span>
868
{| class="formulaSCP" style="width: 100%; text-align: left;" 
869
|-
870
| 
871
{| style="text-align: left; margin:auto;" 
872
|-
873
| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \Vert {}^{n+1}\bar {v}^{i+1}- {}^{n+1}\bar{v}^{i}\Vert \le e_v \Vert {}^{n}\bar {v}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {p}^{i+1}-{}^{n+1}\bar {p}^i\Vert \le e_p  \Vert {}^{n}\bar{p}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {T}^{i+1} - {}^{n+1}\bar {T}^i\Vert \le e_T  \Vert {}^{n}\bar{T}\Vert  \end{array}   </math>
874
|}
875
| style="width: 5px;text-align: right;" | (55)
876
|}
877
878
where <math display="inline">e_v</math>, <math display="inline">e_p</math> and <math display="inline">e_T</math> are prescribed error norms. In the examples presented in the paper we have set <math display="inline">e_v = e_p=e_T = 10^{-3}</math>.
879
880
If  conditions ([[#eq-55|55]]) are satisfied then make <math display="inline">n \leftarrow n+1 </math> and proceed to the next  time step. Otherwise, make the iteration counter <math display="inline">i \leftarrow i+1 </math> and repeat Steps 1&#8211;9.
881
882
The derivatives and integrals in all the matrices of the iteration matrix <math display="inline">{H}_v</math>  are computed on the ''updated discretized configuration at time'' <math display="inline">n</math> (<math display="inline">{}^nV</math>) while the residual vectors <math display="inline">\bar{r}_T</math> and <math display="inline">\bar{r}_m</math> are computed in the current configuration. This is a particular version of the ''updated Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-10'></span><span id='citeF-37'></span><span id='citeF-47'></span><span id='citeF-48'></span>[[#cite-3|[3,10,37,47,48]]].
883
884
'''Remark 7.''' The  tangent deviatoric constitutive matrix <math display="inline">{\boldsymbol{\cal C}}_T</math> in Eq.([[#eq-42|42]]) depends on the constitutive model chosen. For Newtonian fluids <math display="inline">{\boldsymbol{\cal C}}_T = {\boldsymbol{\cal C}}</math> where <math display="inline">{\boldsymbol{\cal C}}</math> is obtained by transforming the deviatoric constitutive tensor <math display="inline">{c}_{ijkl}=\mu (\delta _{ik}\delta _{jl} + \delta _{il}\delta _{jk})</math> from the updated to the current configuration as <math display="inline"> {\cal C}_{ijkl}=F^{-1}_{A_i}F^{-1}_{B_j} F^{-1}_{C_k} F^{-1}_{D_l}     {c}_{ABCD}</math> <span id='citeF-3'></span><span id='citeF-37'></span>[[#cite-3|[3,37]]]. Matrix <math display="inline">\boldsymbol{\cal C}</math> is formed by applying Voigt rule to the terms of tensor <math display="inline">{\cal C}_{ijkl}</math>. A similar expression is obtained for elastic solids changing <math display="inline">\mu </math> by <math display="inline">\bar \mu \Delta t</math> <span id='citeF-10'></span>[[#cite-10|[10]]].
885
886
The second term in the integral of Eq.([[#eq-42|42]]) represents the initial stress contribution to matrix <math display="inline">{K}</math> <span id='citeF-3'></span><span id='citeF-10'></span><span id='citeF-37'></span>[[#cite-3|[3,10,37]]]. Matrices <math display="inline">{G}</math> and <math display="inline">[ S]</math> in Eq.([[#eq-42|42]]) are given by (for 2D problems)
887
888
{| class="formulaSCP" style="width: 100%; text-align: left;" 
889
|-
890
| 
891
{| style="text-align: left; margin:auto;" 
892
|-
893
| style="text-align: center;" | <math> {G} = \begin{bmatrix}     \bar {G} & \bar {0}\\ \bar {0}& \bar {G}     \end{bmatrix}~~,~~  \bar {G} = \begin{bmatrix}   {}_nN_{1,1} & 0 &|& {}_nN_{2,1} & 0 &|& \cdots & {}_nN_{n,1}\\     {}_nN_{1,2} & 0 &|& {}_nN_{2,2} & 0 &|& \cdots & {}_nN_{n,1}     \end{bmatrix} ~~,~~  \bar {0} = \left\{\begin{matrix} 0\\0\end{matrix} \right\}     </math>
894
|}
895
|}
896
897
{| class="formulaSCP" style="width: 100%; text-align: left;" 
898
|-
899
| 
900
{| style="text-align: left; margin:auto;" 
901
|-
902
| style="text-align: center;" | <math>[{ S}] =  \begin{bmatrix}{S} &  {0}\\ {0}& {S}     \end{bmatrix}~~,~~ {0} = \begin{bmatrix}0&0\\0&0\end{bmatrix} ~~,~~     {}_nN_{i,j}={\partial N_i \over \partial ^{n}x_j}     </math>
903
|}
904
| style="width: 5px;text-align: right;" | (56)
905
|}
906
907
where <math display="inline">{S}</math> is the 2d Piola-Kirchhoff stress tensor.
908
909
'''Remark 8.''' The iteration matrix <math display="inline">\hbox{H}_v</math> in Eq.([[#eq-40|40]]) is an ''approximation of the exact tangent matrix for the solid/fluid problem'' in the updated Lagrangian formulation <span id='citeF-3'></span><span id='citeF-10'></span><span id='citeF-37'></span>[[#cite-3|[3,10,37]]]. The form of <math display="inline">\hbox{H}_v</math>  used in this work has yielded good results  with convergence achieved for the nodal values for the velocities,  the pressure and the temperature in 3&#8211;4 iterations for all the problems analyzed.
910
911
'''Remark 9.'''  Including the  bulk stiffness matrix <math display="inline">\hbox{K}_v</math> in <math display="inline">\hbox{H}_v</math> is important for the fast convergence, mass preservation and overall accuracy of the iterative solution for quasi and fully incompressible materials <span id='citeF-9'></span><span id='citeF-38'></span>[[#cite-9|[9,38]]].     The element expression of <math display="inline">\hbox{K}_v</math> can be obtained as <span id='citeF-38'></span>[[#cite-38|[38]]]
912
913
<span id="eq-57"></span>
914
{| class="formulaSCP" style="width: 100%; text-align: left;" 
915
|-
916
| 
917
{| style="text-align: left; margin:auto;" 
918
|-
919
| style="text-align: center;" | <math>\hbox{K}_v^e = \int _{{}^nV^e} \hbox{B}^T \hbox{m} \theta \Delta t\kappa \hbox{m}^T \hbox{B}dV </math>
920
|}
921
| style="width: 5px;text-align: right;" | (57)
922
|}
923
924
where <math display="inline">\theta </math> is a positive number such that <math display="inline">0<   \theta \le 1</math> that has the role of preventing the ill-conditioning of the iteration matrix <math display="inline">{H}_v</math> for highly incompressible materials. For a fully incompressible material (<math display="inline">\kappa =\infty </math>), a finite  value of <math display="inline">\kappa </math> is used in practice within <math display="inline">{K}_v</math> as this helps to obtaining an accurate solution with reduced mass loss in few iterations per time step <span id='citeF-9'></span>[[#cite-9|[9]]]. These considerations, however, do not affect the value of <math display="inline">\kappa </math> within matrix <math display="inline">{M}_1</math> in Eq.([[#eq-36|36]]) that vanishes for a fully incompressible material. A similar approach for improving mass conservation in incompressible fluids was proposed  in <span id='citeF-41'></span>[[#cite-41|[41]]].
925
926
'''Remark 10.''' The time step within a time interval <math display="inline">[n,n+1]</math> is chosen as <math display="inline">\Delta t =\min \left(\frac{^n l_{\min }^e}{\vert {}^n{v}\vert _{\max }},\Delta t_b\right)</math> where <math display="inline">^n l_{\min }^e</math> is the minimum characteristic distance  of all elements in the mesh, with <math display="inline">l^e</math> computed as explained in Section 3, <math display="inline">\vert ^n{v}\vert _{\max }</math> is the maximum value of the modulus of the velocity of all nodes in the mesh  and <math display="inline">\Delta t_b</math> is the critical time step of all nodes approaching a solid boundary defined as <math display="inline">\Delta t_b = \min \left(\frac{^n l_b}{\vert ^n{v}_b\vert _{\max }}\right)</math> where <math display="inline">^n l_b</math> is the distance from the node to the boundary and <math display="inline">^n{v}_b</math> is the velocity of the node. This definition of <math display="inline">\Delta t</math> intends that no node crosses a solid boundary during a time step <span id='citeF-38'></span>[[#cite-38|[38]]].
927
928
'''Remark 11.''' The material properties for the fluid and the solid may be dependent on the temperature. This effect is accounted for by updating the material properties in terms of the temperature within the iteration loop.
929
930
==7 ABOUT THE PARTICLE FINITE ELEMENT METHOD (PFEM)==
931
932
===7.1 The basis of the PFEM===
933
934
Let us consider a domain <math display="inline">V</math> containing  fluid and solid subdomains. Each subdomain is characterized by a set of points, hereafter termed ''virtual particles''. The virtual particles contain all the information  for defining the geometry and the material and mechanical properties of the underlying subdomain. In the PFEM both subdomains are modelled using an ''updated'' ''Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-47'></span>[[#cite-3|[3,47]]].
935
936
The solution steps within a time step in the PFEM are as follows:
937
938
<div id='img-2'></div>
939
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
940
|-
941
|[[Image:draft_Samper_718004671-Figure2a.png|600px|Sequence of steps to update a “cloud” of virtual particles (nodes) representing a domain containing a fluid and a solid  from ⁿt  to ⁿ⁺²t. u, v,a,ɛ , ̇ɛ and σ  denote the displacement, the velocity, the acceleration, the strain, the strain rates and the Cauchy stresses, respectively.]]
942
|- style="text-align: center; font-size: 75%;"
943
| colspan="1" | '''Figure 2:''' Sequence of steps to update a “cloud” of virtual particles (nodes) representing a domain containing a fluid and a solid  from <math>^n t</math>  to <math>{}^{n+2}t</math>. <math>{u}, {v},{a},{\boldsymbol \varepsilon } , \dot {\boldsymbol \varepsilon }</math> and <math>{\boldsymbol \sigma }</math>  denote the displacement, the velocity, the acceleration, the strain, the strain rates and the Cauchy stresses, respectively.
944
|}
945
946
<ol>
947
948
<li>The starting point at each time step is the cloud of points <math display="inline">C</math> in the fluid and solid domains. For instance <math display="inline">^{n} C</math> denotes the cloud at time <math display="inline">t={}^n t </math> (Figure [[#img-1|1]]). The virtual particles are assumed to be coincident with the nodes in the finite element mesh generated in Step 3. </li>
949
950
<li>Identify the boundaries defining the analysis domain <math display="inline">^{n} V</math>, as well as the subdomains in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution, including separation and re-entering of nodes. The Alpha Shape method <span id='citeF-8'></span>[[#cite-8|[8]]] is used for the boundary definition. Clearly, the accuracy in the reconstruction of the boundaries depends on the number of points in the vicinity of each boundary and on the Alpha Shape parameter. In the problems solved in this work the Alpha Shape method has been implementation as described in <span id='citeF-14'></span><span id='citeF-26'></span>[[#cite-14|[14,26]]]. </li>
951
952
<li> Discretize  the analysis domain <math display="inline">^{n} V</math> with a finite element mesh <math display="inline">{}^{n} M</math> using the virtual particles as the mesh nodes. We use an efficient mesh generation scheme based on an enhanced Delaunay tesselation <span id='citeF-13'></span><span id='citeF-14'></span>[[#cite-13|[13,14]]]. </li>
953
954
<li> Solve the Lagrangian equations of motion for the overall continuum using the standard FEM. Compute the state variables in  the next (updated) configuration for <math display="inline">^n t+\Delta t</math>: velocities, pressure, strain rate and viscous stresses in the fluid and displacements, stresses and strains in the solid. </li>
955
956
<li>  Move the mesh nodes to a new position <math display="inline">^{n+1} C</math> where ''n''+1 denotes the time <math display="inline">{}^n t +\Delta t</math>, in terms of the time increment size. </li>
957
958
<li> Go back to step 1 and repeat the solution for the next time step to obtain <math display="inline">{}^{n+2} C</math>. </li>
959
960
</ol>
961
962
Note that the key differences between the PFEM and the classical FEM are the remeshing technique and the identification of the domain boundary at each time step.
963
964
The CPU time required for meshing grows linearly with the number of nodes. As a general rule,  meshing consumes for 3D problems around 15% of the total CPU time per time step, while the solution of the equations (with typically 3 iterations per time step) and the system assembly consume approximately 70% and 15% of the CPU time per time step, respectively. These figures refer to analyses in a single processor Pentium IV PC <span id='citeF-35'></span>[[#cite-35|[35]]]. Considerable speed can be gained using parallel computing techniques.
965
966
Application of the PFEM in fluid and solid mechanics and in fluid-structure interaction problems can be found in <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-20'></span><span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-34'></span><span id='citeF-35'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-41'></span><span id='citeF-43'></span>[[#cite-1|[1,2,4,5,8,9,13,14,15,16,17,18,20,26,27,34,35,37,38,41,43]]] as well in www.cimne.com/pfem.
967
968
===7.2 Treatment of contact conditions in the PFEM===
969
970
Known velocities at boundaries in the PFEM are prescribed in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. Surface tractions are applied to the Neumann part of the boundary, as it is usual in the FEM.
971
972
Contact between fluid particles and fixed boundaries is accounted for by the incompressibility condition which  naturally prevents fluid nodes to penetrate into the solid boundaries <span id='citeF-14'></span><span id='citeF-26'></span><span id='citeF-32'></span><span id='citeF-35'></span>[[#cite-14|[14,26,32,35]]].
973
974
The contact between two solid interfaces is treated by introducing a layer of ''contact elements'' between the two interacting solid interfaces (Figure [[#img-3|3]]). This contact layer is ''automatically created during the mesh generation step'' by prescribing a minimum distance <math display="inline">\left(h_{c} \right)</math> between two solid boundaries. If the distance exceeds the minimum value <math display="inline">\left(h_{c} \right)</math> then the generated elements are treated as fluid elements. Otherwise the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacement is introduced  in <span id='citeF-26'></span><span id='citeF-32'></span><span id='citeF-35'></span>[[#cite-26|[26,32,35]]].
975
976
This algorithm allows us to model complex frictional contact conditions between two or more interacting bodies moving in water in  a simple manner. The algorithm has  been used to model frictional contact situations between rigid or elastic solids in structural mechanics applications, such as soil/rock excavation problems <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]]. The frictional contact algorithm described above has been extended by Oliver ''et al.'' <span id='citeF-20'></span>[[#cite-20|[20]]] for analysis of metal cutting and machining problems.
977
978
<div id='img-3'></div>
979
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 65%;max-width: 100%;"
980
|-
981
|[[Image:draft_Samper_718004671-contact-solid-boundaries.png|390px|Layer of contact elements at a soil-solid interface modelled with the PFEM]]
982
|- style="text-align: center; font-size: 75%;"
983
| colspan="1" | '''Figure 3:''' Layer of contact elements at a soil-solid interface modelled with the PFEM
984
|}
985
986
==8 EXAMPLES==
987
988
===8.1 Numerical examples of manufacturing problems===
989
990
Some examples are presented next to show the capabilities of the method for the simulation of industrial forming processes. All the problems presented involve large strains and big changes in the geometry boundaries. The “particle” description of the continuum assumed in the PFEM introduces the capability of adapting the geometry to the computed displacement solution automatically.
991
992
The numerical solution is computed taking in account thermal and mechanical effects with the solution strategy described in Section 6.
993
994
The changes in the geometry are recognized by the PFEM features. As explained in Section 7.1, the finite element mesh is obtained by  a reconnection of the particles and a remeshing of the domain. A refinement of the critical areas at each time step has been performed. The criterion to detect the critical areas for refining the mesh is based on the evaluation of the plastic energy dissipation, but other mechanical variables can be used for this purpose. Details of the application of the PFEM to solid mechanics problems can be found in <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]].
995
996
===8.2 Extrusion of a steel plate===
997
998
Extrusion is one of the techniques used to reduce the section of pieces of metal and form a reduced piece of the same object. This example considers a piece of steel with the dimensions shown in  Figure [[#img-4|4]] pushed towards a rigid wall die. Two different shapes of the wall die have been analysed. Both models are depicted in Figure [[#img-4|4]]. The symmetry of the problem allows us to consider only the upper part of the piece and a 2D simulation to solve the problem. The material is modelled with a thermo-elastoplastic constitutive law with a Huber-Von Misses yield surface. The constitutive model can be found in <span id='citeF-42'></span>[[#cite-42|[42]]]. The physical properties of the material (steel) considered in the example are shown in Box 2.
999
1000
1001
{| class="wikitable" style="text-align: left; margin: 1em auto;"
1002
1003
|-
1004
| colspan='2' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;" | ''MATERIAL PROPERTIES'':
1005
1006
|-
1007
| style="border-left: 2px solid;" | Young  Modulus 
1008
| style="border-right: 2px solid;" | <math> 206.9 \times 10^9  Pa</math>
1009
|-
1010
| style="border-left: 2px solid;" | Poisson Coefficient 
1011
| style="border-right: 2px solid;" | 0.29 
1012
|-
1013
| style="border-left: 2px solid;" | Yield Stress 
1014
| style="border-right: 2px solid;" | <math>450\times 10^6  Pa</math>
1015
|-
1016
| style="border-left: 2px solid;" | Linear Hardening Modulus &nbsp;~&nbsp;
1017
| style="border-right: 2px solid;" | <math>129.24\times 10^6  Pa</math>
1018
|-
1019
| style="border-left: 2px solid;" | Reference Hardening 
1020
| style="border-right: 2px solid;" | <math>450\times 10^6  Pa</math>
1021
|-
1022
| style="border-left: 2px solid;" | Saturation Hardening 
1023
| style="border-right: 2px solid;" | <math>715\times 10^6  Pa</math>
1024
|-
1025
| style="border-left: 2px solid;" | Hardening Exponent 
1026
| style="border-right: 2px solid;" | 16.93 
1027
|-
1028
| style="border-left: 2px solid;" | Conductivity 
1029
| style="border-right: 2px solid;" | <math>45 N/s K</math>
1030
|-
1031
| style="border-left: 2px solid;" | Density 
1032
| style="border-right: 2px solid;" | <math>7800 kg/m^3 </math>
1033
|-
1034
| style="border-left: 2px solid;" | Specific Capacity 
1035
| style="border-right: 2px solid;" | <math>460  m^2/s^2 K</math>
1036
|-
1037
| style="border-left: 2px solid;" | Flow Stress Softening 
1038
| style="border-right: 2px solid;" | <math>0.002 K^{-1}</math>
1039
|-
1040
| style="border-left: 2px solid;" | Hardening Softening 
1041
| style="border-right: 2px solid;" | <math>0.002 K^{-1}</math>
1042
|-
1043
| style="border-left: 2px solid;" | Dissipation Factor 
1044
| style="border-right: 2px solid;" | 0.9 
1045
|- style="border-bottom: 2px solid;"
1046
| style="border-left: 2px solid;" | Expansion Coefficient 
1047
| style="border-right: 2px solid;" | <math>10^{-5} K^{-1} </math>
1048
1049
|}<br/>
1050
1051
'''Box 2.''' Material properties for the steel extrusion problem
1052
1053
An imposed velocity is set on the left side of the models: <math display="inline">10 mm/s</math> for the vertical die wall and <math display="inline">2.5 mm/s</math> for the inclined wall. In both models the vertical displacement is imposed to a zero value at the lower side of the steel plate. The initial temperature is <math display="inline">^\circ </math>T=293.15 K and the left boundary is fixed to this temperature during the analysis.
1054
1055
<div id='img-4'></div>
1056
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1057
|-
1058
|[[Image:draft_Samper_718004671-dimensions_wall_step.png|570px|Extrusion of a steel plate with two wall dies. Dimensions expressed in mm. ]]
1059
|- style="text-align: center; font-size: 75%;"
1060
| colspan="1" | '''Figure 4:''' Extrusion of a steel plate with two wall dies. Dimensions expressed in mm. 
1061
|}
1062
1063
<div id='img-5'></div>
1064
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1065
|-
1066
|[[Image:draft_Samper_718004671-extrusion_step_mesh.png|570px|Extrusion of a steel plate at three  time instants: t = 0.09s, t = 2.59s, t = 5.19s. Vertical die wall]]
1067
|- style="text-align: center; font-size: 75%;"
1068
| colspan="1" | '''Figure 5:''' Extrusion of a steel plate at three  time instants: t = 0.09s, t = 2.59s, t = 5.19s. Vertical die wall
1069
|}
1070
1071
The results of the extrusion analysis with the vertical step are shown in Figures [[#img-5|5]] and [[#img-6|6]]. Different time instants are considered, from the beginning to the end of the process. The contact forces with the wall, the temperature of the model, the plastic strain and the plastic strain rate are measures that can be obtained from the numerical results. With the distribution and the values of these quantities the extrusion process can be studied and optimized.
1072
1073
A similar set of results for a different extrusion problems using a die wall with a progressive reduction of the section are depicted in Figures [[#img-7|7]] and [[#img-8|8]].
1074
1075
<div id='img-6'></div>
1076
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1077
|-
1078
|[[Image:draft_Samper_718004671-extrusion_step_results.png|570px|Results of the extrusion of a steel plate  (vertical die wall) at t = 2.59s. Nodal forces in Newtons, pressure in N/mm², temperature in Kelvins and plastic strain rate in 1/s. ]]
1079
|- style="text-align: center; font-size: 75%;"
1080
| colspan="1" | '''Figure 6:''' Results of the extrusion of a steel plate  (vertical die wall) at t = 2.59s. Nodal forces in Newtons, pressure in <math>N/mm^2</math>, temperature in Kelvins and plastic strain rate in <math>1/s</math>. 
1081
|}
1082
1083
<div id='img-7'></div>
1084
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1085
|-
1086
|[[Image:draft_Samper_718004671-extrusion_slope_mesh.png|570px|Extrusion of a steel plate at three  time instants: t = 0.9s, t = 11.9s, t = 24.9s. Inclined die wall]]
1087
|- style="text-align: center; font-size: 75%;"
1088
| colspan="1" | '''Figure 7:''' Extrusion of a steel plate at three  time instants: t = 0.9s, t = 11.9s, t = 24.9s. Inclined die wall
1089
|}
1090
1091
<div id='img-8'></div>
1092
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1093
|-
1094
|[[Image:draft_Samper_718004671-extrusion_slope_results.png|570px|Results of the extrusion of a steel plate at t = 11.9s (inclined die wall). Nodal forces in Newtons, pressure in N/mm², temperature in Kelvins and plastic strain rate in 1/s. ]]
1095
|- style="text-align: center; font-size: 75%;"
1096
| colspan="1" | '''Figure 8:''' Results of the extrusion of a steel plate at t = 11.9s (inclined die wall). Nodal forces in Newtons, pressure in <math>N/mm^2</math>, temperature in Kelvins and plastic strain rate in <math>1/s</math>. 
1097
|}
1098
1099
<div id='img-9'></div>
1100
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1101
|-
1102
|[[Image:draft_Samper_718004671-single_shear_cut.png|570px|Cutting of a steel plate with a single cutter. Dimensions in mm.]]
1103
|- style="text-align: center; font-size: 75%;"
1104
| colspan="1" | '''Figure 9:''' Cutting of a steel plate with a single cutter. Dimensions in mm.
1105
|}
1106
1107
<div id='img-10'></div>
1108
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1109
|-
1110
|[[Image:draft_Samper_718004671-double_shear_cut.png|570px|Cutting of a steel plate with a two symmetric cutters. Dimensions in mm.]]
1111
|- style="text-align: center; font-size: 75%;"
1112
| colspan="1" | '''Figure 10:''' Cutting of a steel plate with a two symmetric cutters. Dimensions in mm.
1113
|}
1114
1115
===8.3 Shear cutting===
1116
1117
The next example shows the cutting of solid plate in 2D. It is a thermo-mechanical problem which involves large deformations and the process of cutting by means of a rigid tool. The initial conditions of the problem are depicted in Figures [[#img-9|9]] and [[#img-10|10]]. A thermo-elastoplastic constitutive law with the same  material properties as described in the previous example has been used. As previously mentioned, the PFEM combines the Lagrangian updating of the configuration with a redefinition of the domain geometry. These features allow to model naturally the cutting region by  refining  the solid model in the tool tip and by  inserting  new particles (nodes) when large flat boundaries appear at the tip contact zone.
1118
1119
The cutting tool moves downwards with a velocity of <math display="inline">5 mm/s</math> in the single cutter model and with a velocity of <math display="inline">200mm/s</math> in the problem with two cutters. In both models the horizontal displacement is imposed to zero in the outer sides of the plate. The initial temperature in the domain is set to 293.15 K and is kept fixed to this value in the outer boundary sides during  the analysis.
1120
1121
The PFEM results of the material deformation and the temperature distribution for various time instants are shown in Figures [[#img-11|11]] and [[#img-12|12]].
1122
1123
<div id='img-11'></div>
1124
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1125
|-
1126
|[[Image:draft_Samper_718004671-single_cut_results.png|570px|Cutting process. Results at the onset,  the middle and  the end of the process using a single cutter. Nodal forces are expressed in Newtons and the temperature in Kelvins.]]
1127
|- style="text-align: center; font-size: 75%;"
1128
| colspan="1" | '''Figure 11:''' Cutting process. Results at the onset,  the middle and  the end of the process using a single cutter. Nodal forces are expressed in Newtons and the temperature in Kelvins.
1129
|}
1130
1131
<div id='img-12'></div>
1132
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1133
|-
1134
|[[Image:draft_Samper_718004671-double_cut_results.png|570px|Results at the onset, the middle and  the end of a plate cutting process with two symmetric cutters. Nodal forces are expressed in Newtons and the temperature in Kelvins.]]
1135
|- style="text-align: center; font-size: 75%;"
1136
| colspan="1" | '''Figure 12:''' Results at the onset, the middle and  the end of a plate cutting process with two symmetric cutters. Nodal forces are expressed in Newtons and the temperature in Kelvins.
1137
|}
1138
1139
<div id='img-13'></div>
1140
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1141
|-
1142
|[[Image:draft_Samper_718004671-dimensions_forge.png|570px|Forging of a steel cylinder against to curved rigid walls. Dimensions in meters.]]
1143
|- style="text-align: center; font-size: 75%;"
1144
| colspan="1" | '''Figure 13:''' Forging of a steel cylinder against to curved rigid walls. Dimensions in meters.
1145
|}
1146
1147
===8.4 Forging of a cylindrical steel piece===
1148
1149
The next example considers the forging of a  steel cylinder. The material properties are the same as in Box 2. The mould is defined by two rigid walls that move  towards each other at the same speed, (<math display="inline">1 m/s</math>)  and enclose the steel piece. The dimensions of the model are depicted in Figure [[#img-13|13]]. An axisymmetric solution is considered. The initial temperature of the steel material is <math display="inline">^\circ </math>T =293.15 K.
1150
1151
The results of the material deformation, the pressure contours and the temperature distribution for various time instants are shown in Figures [[#img-14|14]] and [[#img-15|15]].
1152
1153
<div id='img-14'></div>
1154
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1155
|-
1156
|[[Image:draft_Samper_718004671-forge_mesh_force.png|570px|Meshes and wall nodal forces at the beginning,  the middle and  the end of the forging process. Nodal forces are expressed in Newtons.]]
1157
|- style="text-align: center; font-size: 75%;"
1158
| colspan="1" | '''Figure 14:''' Meshes and wall nodal forces at the beginning,  the middle and  the end of the forging process. Nodal forces are expressed in Newtons.
1159
|}
1160
1161
<div id='img-15'></div>
1162
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 95%;max-width: 100%;"
1163
|-
1164
|[[Image:draft_Samper_718004671-forge_pressure_temperature.png|570px|Pressure and temperatures at the beginning,  the middle and  the end of the forging process. Pressure is expressed in Pascals and the temperature in Kelvins.]]
1165
|- style="text-align: center; font-size: 75%;"
1166
| colspan="1" | '''Figure 15:''' Pressure and temperatures at the beginning,  the middle and  the end of the forging process. Pressure is expressed in Pascals and the temperature in Kelvins.
1167
|}
1168
1169
===8.5 Hot forging of a metal piece===
1170
1171
Metal hot forging consists on deforming a work piece  at high temperature under compressive forces in order to obtain the desired shape of the final product. It is a largely used process in industrial metal manufacturing as it allows an improvement of the material properties of the final part via a controlled deformation process. For example, via  hot forging it is possible to reduce the impurities and  defects in the material, thereby avoiding stress concentrations in the final product.
1172
1173
The metal forging can be performed using an open die or a closed die. In this work, the latter case is considered. The geometry and the problem data for the 2D simulation are shown in Figure [[#img-16|16]]. The forging process is performed by controlling the displacement of the upper rectangular die which is pushed downwards with a velocity v=0.01<math display="inline">m/s</math> for a duration of 5.65<math display="inline">s</math>, while the material rigid walls are kept at the same position during the analysis. At the beginning of the process the metal piece and the lower walls have temperatures <math display="inline">^\circ </math>T=2000 K and <math display="inline">^\circ </math>T =1000 K, respectively. The forging is carried out through two steps: during the first one the piece is compressed by the descending rectangular  die; then the die is stopped and from t=6.0<math display="inline">s</math> to t=50.0<math display="inline">s</math> the piece is cooled down.  The cooling process is carried out by applying the temperature given by the following step function to both the upper rectangular die and the internal walls.
1174
1175
<span id="eq-58"></span>
1176
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1177
|-
1178
| 
1179
{| style="text-align: left; margin:auto;" 
1180
|-
1181
| style="text-align: center;" | <math>\displaystyle  \begin{array}{l}\displaystyle  T=850K   6.0s < t \leq 8.5s\\ \displaystyle T=600K   8.5s < t \leq 10.5s\\ \displaystyle T=300K   10.5s < t \leq 50s\\ \end{array}   </math>
1182
|}
1183
| style="width: 5px;text-align: right;" | (58)
1184
|}
1185
1186
It is assumed that the rectangular die has the same temperature of the metal piece during the compression phase. The normal heat flux through the surfaces in contact with the air has not been taken into account in the analysis.
1187
1188
The metal is modelled as an isotropic incompressible non-Newtonian fluid with a non-linear viscosity  given by the following relation <span id='citeF-44'></span><span id='citeF-45'></span>[[#cite-44|[44,45]]]
1189
1190
<span id="eq-59"></span>
1191
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1192
|-
1193
| 
1194
{| style="text-align: left; margin:auto;" 
1195
|-
1196
| style="text-align: center;" | <math>\mu = \frac{\sigma _{y}}{\sqrt{3}\bar \varepsilon }  </math>
1197
|}
1198
| style="width: 5px;text-align: right;" | (59)
1199
|}
1200
1201
where <math display="inline">{\sigma _{y}}</math> is the uniaxial yield stress of the material and <math display="inline">\bar{{\varepsilon }}</math> is the deviatoric strain rate invariant defined as:
1202
1203
<span id="eq-60"></span>
1204
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1205
|-
1206
| 
1207
{| style="text-align: left; margin:auto;" 
1208
|-
1209
| style="text-align: center;" | <math>\bar \varepsilon = \sqrt{\frac{2}{3}{\varepsilon }_{ij} {\varepsilon }_{ij}}  </math>
1210
|}
1211
| style="width: 5px;text-align: right;" | (60)
1212
|}
1213
1214
The uniaxial yield stress <math display="inline">{\sigma _{y}}</math> depends on the temperature as follows:
1215
1216
<span id="eq-61"></span>
1217
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1218
|-
1219
| 
1220
{| style="text-align: left; margin:auto;" 
1221
|-
1222
| style="text-align: center;" | <math>\displaystyle  \begin{array}{l}\displaystyle {\sigma _{y}}=157.7                                                            T \leq 1500\\ \displaystyle {\sigma _{y}}=157.7+3.6\cdot (T-1500)^2-3250\cdot (T-1500)     1500< T \leq 1930 \\ \displaystyle {\sigma _{y}}=79.4                                                               T >1930 \\ \end{array}   </math>
1223
|}
1224
| style="width: 5px;text-align: right;" | (61)
1225
|}
1226
1227
The constitutive equation ([[#eq-59|59]]) is typical in the study of industrial forming processes by the plastic/viscoplastic flow approach <span id='citeF-44'></span><span id='citeF-45'></span><span id='citeF-48'></span>[[#cite-44|[44,45,48]]].
1228
1229
The metal piece has been discretized using 6849 3-noded triangular finite elements. The 2D coupled thermal-mechanical simulation has been run for 50<math display="inline">s</math> with <math display="inline">\Delta t</math>=<math display="inline">0.0005s</math>.
1230
1231
The snapshots of Figure [[#img-17|17]] refer to the initial compression phase that has a duration of 5.65<math display="inline">s</math>. The upper rectangular die moves downwards forcing the metal piece to take the shape of the lower rigid walls.
1232
1233
Figure [[#img-18|18]] shows four representative snapshots of the cooling period.
1234
1235
<div id='img-16'></div>
1236
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1237
|-
1238
|[[Image:draft_Samper_718004671-forgingInput2.png|400px|]]
1239
|[[Image:draft_Samper_718004671-forgingData.png|400px|]]
1240
|-
1241
| colspan="2"|[[Image:draft_Samper_718004671-metalData.png|600px|Hot forging of a metal piece. Initial geometry, coordinates of corner points and material properties.]]
1242
|- style="text-align: center; font-size: 75%;"
1243
| colspan="2" | '''Figure 16:''' Hot forging of a metal piece. Initial geometry, coordinates of corner points and material properties.
1244
|}
1245
1246
<div id='img-17'></div>
1247
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1248
|-
1249
|[[Image:draft_Samper_718004671-forging800.png|400px|]]
1250
|[[Image:draft_Samper_718004671-forging1500.png|400px|]]
1251
|-
1252
|[[Image:draft_Samper_718004671-forging3610.png|400px|]]
1253
|[[Image:draft_Samper_718004671-tempComp.png|400px|]]
1254
|-
1255
|[[Image:draft_Samper_718004671-forging5320.png|400px|]]
1256
|[[Image:draft_Samper_718004671-forging5560.png|400px|]]
1257
|-
1258
| colspan="2"|[[Image:draft_Samper_718004671-forging5650.png|300px|Hot forging of a metal piece. Compression phase. Snapshots of the deformation with temperature contours at different time instants.]]
1259
|- style="text-align: center; font-size: 75%;"
1260
| colspan="2" | '''Figure 17:''' Hot forging of a metal piece. Compression phase. Snapshots of the deformation with temperature contours at different time instants.
1261
|}
1262
1263
<div id='img-18'></div>
1264
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1265
|-
1266
|[[Image:draft_Samper_718004671-forgingCooling8450.png|400px|]]
1267
|[[Image:draft_Samper_718004671-forgingCooling10445.png|400px|]]
1268
|-
1269
|[[Image:draft_Samper_718004671-forgingCooling25000.png|400px|]]
1270
|[[Image:draft_Samper_718004671-forgingCooling50000.png|400px|]]
1271
|-
1272
| colspan="2"|[[Image:draft_Samper_718004671-tempCooling.png|300px|Hot forging of a metal piece. Cooling period. Snapshots of the  temperature contours on the deformed geometry at different time steps.]]
1273
|- style="text-align: center; font-size: 75%;"
1274
| colspan="2" | '''Figure 18:''' Hot forging of a metal piece. Cooling period. Snapshots of the  temperature contours on the deformed geometry at different time steps.
1275
|}
1276
1277
<div id='img-19'></div>
1278
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1279
|-
1280
|[[Image:draft_Samper_718004671-ballsInput2.png|400px|]]
1281
|[[Image:draft_Samper_718004671-ballDataGeo.png|400px|]]
1282
|-
1283
| colspan="2"|[[Image:draft_Samper_718004671-ballDataMat.png|600px|Falling of three solid objects in a heated tank filled with a fluid. Initial geometry, thermal conditions and material properties.]]
1284
|- style="text-align: center; font-size: 75%;"
1285
| colspan="2" | '''Figure 19:''' Falling of three solid objects in a heated tank filled with a fluid. Initial geometry, thermal conditions and material properties.
1286
|}
1287
1288
===8.6 Falling of three solid objects in a heated tank filled with fluid===
1289
1290
Three solid objects with the same shape fall from the same height into a tank containing a fluid at rest. The geometry, the problem data, and the initial thermal conditions are shown in Figure [[#img-19|19]].
1291
1292
The fluid in the tank has an initial temperature of T=340<math display="inline">K</math>, while the solid bodies from left to right have initial temperatures of T=180 K, T=200 K, and T=220 K, respectively. The solid and the fluid domains have been discretized with a finite element mesh of 9394 3-noded triangular elements. The simulation has been run for eight seconds with <math display="inline">\Delta t</math>=<math display="inline">0.0001 s</math>. The heat flux in the normal direction is assumed to be zero for the boundaries in contact with the air and the walls. Figure [[#img-20|20]] shows six representative snapshots of the numerical simulation with the temperature results plotted over the fluid domain and the objects domains.
1293
1294
<div id='img-20'></div>
1295
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1296
|-
1297
|[[Image:draft_Samper_718004671-Balls664.png|400px|]]
1298
|[[Image:draft_Samper_718004671-Balls1816.png|400px|]]
1299
|-
1300
|[[Image:draft_Samper_718004671-Balls2656.png|400px|]]
1301
|[[Image:draft_Samper_718004671-tempIce.png|400px|]]
1302
|-
1303
|[[Image:draft_Samper_718004671-Balls4500.png|400px|]]
1304
|[[Image:draft_Samper_718004671-Balls7000.png|400px|]]
1305
|-
1306
| colspan="2"|[[Image:draft_Samper_718004671-Balls8000.png|300px|Falling of three solid objects in a heated tank filled with a fluid. Snapshots with temperature contours at different time steps.]]
1307
|- style="text-align: center; font-size: 75%;"
1308
| colspan="2" | '''Figure 20:''' Falling of three solid objects in a heated tank filled with a fluid. Snapshots with temperature contours at different time steps.
1309
|}
1310
1311
<div id='img-21'></div>
1312
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 85%;max-width: 100%;"
1313
|-
1314
|[[Image:draft_Samper_718004671-REPBalls6hGraph.png|510px|Falling of three solid objects in a heated tank filled with a fluid. Evolution of the temperature at the center of the three objects.]]
1315
|- style="text-align: center; font-size: 75%;"
1316
| colspan="1" | '''Figure 21:''' Falling of three solid objects in a heated tank filled with a fluid. Evolution of the temperature at the center of the three objects.
1317
|}
1318
1319
In the graph of Figure [[#img-21|21]] the evolution of the temperature at the central point of the three solid objects is plotted.
1320
1321
===8.7 Melting of an ice block===
1322
1323
A cylindrical ice block at initial temperature <math display="inline">^\circ </math>T=270 K is dropped into a tank containing water at rest at an initial temperature <math display="inline">^\circ </math>T=340 K. The initial geometry, the problem data and the initial thermal conditions for the 2D simulation are shown in Figure [[#img-22|22]].
1324
1325
<div id='img-22'></div>
1326
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1327
|-
1328
|[[Image:draft_Samper_718004671-iceInput2.png|400px|]]
1329
|[[Image:draft_Samper_718004671-iceDataGeo.png|400px|]]
1330
|-
1331
| colspan="2"|[[Image:draft_Samper_718004671-iceDataMat.png|600px|Melting of an ice cylinder. Geometry, material data and initial thermal conditions.]]
1332
|- style="text-align: center; font-size: 75%;"
1333
| colspan="2" | '''Figure 22:''' Melting of an ice cylinder. Geometry, material data and initial thermal conditions.
1334
|}
1335
1336
Ice is treated as a hypoelastic solid until  some of its elements reach the fusion temperature (T=273.15 K ). These elements are transferred to the fluid domain and take the physical properties of the fluid. For this analysis the following assumptions were made: the mechanical and thermal properties of water and  ice do not change with  temperature and the heat normal flux along the boundaries in contact with the air or the walls have been considered to be zero.
1337
1338
In Figure [[#img-23|23]] snapshots of some representative instants of the analysis are shown. The temperature contours are plotted over the water domain and the ice block.
1339
1340
In Figure [[#img-24|24]] the detail of the melting of the ice block is illustrated. The finite element mesh is drawn over the solid and the fluid domains.
1341
1342
<div id='img-23'></div>
1343
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1344
|-
1345
|[[Image:draft_Samper_718004671-Ice44.png|400px|]]
1346
|[[Image:draft_Samper_718004671-Ice110.png|400px|]]
1347
|-
1348
|[[Image:draft_Samper_718004671-Ice263.png|400px|]]
1349
|[[Image:draft_Samper_718004671-tempIce.png|400px|]]
1350
|-
1351
|[[Image:draft_Samper_718004671-Ice455.png|400px|]]
1352
|[[Image:draft_Samper_718004671-Ice1001.png|400px|]]
1353
|-
1354
| colspan="2"|[[Image:draft_Samper_718004671-Ice1600.png|300px|Melting of an ice cylinder. Snapshots of the melting process with temperature contours at different time steps.]]
1355
|- style="text-align: center; font-size: 75%;"
1356
| colspan="2" | '''Figure 23:''' Melting of an ice cylinder. Snapshots of the melting process with temperature contours at different time steps.
1357
|}
1358
1359
<div id='img-24'></div>
1360
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1361
|-
1362
|[[Image:draft_Samper_718004671-iceZoom281.png|400px|]]
1363
|[[Image:draft_Samper_718004671-iceZoom470.png|400px|]]
1364
|-
1365
|[[Image:draft_Samper_718004671-iceZoom977.png|400px|]]
1366
|[[Image:draft_Samper_718004671-iceZoom1274.png|400px|]]
1367
|-
1368
|[[Image:draft_Samper_718004671-iceZoom1469.png|400px|]]
1369
|[[Image:draft_Samper_718004671-iceZoom1547.png|400px|Melting of an ice cylinder. Zoom of the ice domain at different time instants.]]
1370
|- style="text-align: center; font-size: 75%;"
1371
| colspan="2" | '''Figure 24:''' Melting of an ice cylinder. Zoom of the ice domain at different time instants.
1372
|}
1373
1374
===8.8 Falling and warming of a solid sphere in a cilindrical tank containing hot water===
1375
1376
A solid sphere at <math display="inline">^\circ </math>T=270 K is dropped into a cilndrical tank containing still water  at <math display="inline">^\circ </math>T=340 K. In  Figure [[#img-22|22]] the geometry and the problem data are given. The 3D simulation is run for a duration of 3<math display="inline">s</math>  with <math display="inline">\Delta t</math>=<math display="inline">0.001 s</math>. A finite element mesh of 144663 4-noded tetrahedral elements has been used to discretize both the fluid domain and the sphere. Figure [[#img-26|26]] shows six representative snapshots of the analysis. The temperature contours have been plotted over the fluid domain and the sphere.
1377
1378
Figure [[#img-27|27]] shows the  evolution of the temperature at the center of the sphere versus time. The solid material has high conductivity and this explains its fast warming. After three seconds, it practically reaches the equilibrium temperature. We note that the  normal heat flux along the boundaries in contact with the air have been considered to be zero.
1379
1380
<div id='img-25'></div>
1381
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1382
|-
1383
|[[Image:draft_Samper_718004671-ball3DInput.png|400px|]]
1384
|[[Image:draft_Samper_718004671-ball3DDataGeo.png|400px|]]
1385
|-
1386
| colspan="2"|[[Image:draft_Samper_718004671-ball3DDataMat.png|600px|Warming of a solid sphere with initial temperature ^∘T=270 K falling into a cylindrical tank containing a fluid at ^∘T=340 K. Initial geometry and material properties.]]
1387
|- style="text-align: center; font-size: 75%;"
1388
| colspan="2" | '''Figure 25:''' Warming of a solid sphere with initial temperature <math>^\circ </math>T=270 K falling into a cylindrical tank containing a fluid at <math>^\circ </math>T=340 K. Initial geometry and material properties.
1389
|}
1390
1391
<div id='img-26'></div>
1392
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1393
|-
1394
|[[Image:draft_Samper_718004671-ball3D8.png|400px|]]
1395
|[[Image:draft_Samper_718004671-ball3D17.png|400px|]]
1396
|-
1397
|[[Image:draft_Samper_718004671-ball3D38.png|400px|]]
1398
|[[Image:draft_Samper_718004671-tempIce.png|400px|]]
1399
|-
1400
|[[Image:draft_Samper_718004671-ball3D53.png|400px|]]
1401
|[[Image:draft_Samper_718004671-ball3D71.png|400px|]]
1402
|-
1403
| colspan="2"|[[Image:draft_Samper_718004671-ball3D300.png|300px|Warming of a solid sphere (^∘T=270 K) falling into a cylindrical tank containing a fluid at ^∘T=340 K. Snapshots of the  motion of the sphere and the temperature contours at different time steps.]]
1404
|- style="text-align: center; font-size: 75%;"
1405
| colspan="2" | '''Figure 26:''' Warming of a solid sphere (<math>^\circ </math>T=270 K) falling into a cylindrical tank containing a fluid at <math>^\circ </math>T=340 K. Snapshots of the  motion of the sphere and the temperature contours at different time steps.
1406
|}
1407
1408
<div id='img-27'></div>
1409
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
1410
|-
1411
|[[Image:draft_Samper_718004671-ball3DEvolution.png|480px|Warming of a solid sphere (^∘T=270 K) falling into a cylindrical tank containing a fluid at ^∘T=340 K. Evolution of the temperature at the center of the sphere.]]
1412
|- style="text-align: center; font-size: 75%;"
1413
| colspan="1" | '''Figure 27:''' Warming of a solid sphere (<math>^\circ </math>T=270 K) falling into a cylindrical tank containing a fluid at <math>^\circ </math>T=340 K. Evolution of the temperature at the center of the sphere.
1414
|}
1415
1416
==9 CONCLUDING REMARKS==
1417
1418
We have presented a unified Lagrangian formulation for analysis of industrial forming processes involving thermally coupled interactions between deformable continua containing fluids and solids. A  residual-based  expression of the mass conservation equation obtained using the FIC method  provides the necessary stability for quasi/fully incompressible situations. The  governing equations for the generalized continuum are discretized with the FEM using simplicial elements with equal linear interpolation for the velocities, the pressure and the temperature. The merits of the formulation in terms of its general applicability have been demonstrated in the solution of a variety of thermally-coupled industrial forming processes using the Particle Finite Element Method (PFEM).
1419
1420
==ACKNOWLEDGEMENTS==
1421
1422
This research was partially supported by the Advanced Grant  project SAFECON of the European Research Council and the HFLUIDS project of the National Research Programme of Spain.
1423
1424
===BIBLIOGRAPHY===
1425
1426
<div id="cite-1"></div>
1427
'''[[#citeF-1|[1]]]'''  Aubry R, Idelsohn SR, Oñate E (2005)  Particle finite element method in fluid-mechanics including thermal convection-diffusion. Computers and Structures,  Vol. 83, (17-18), pp 1459&#8211;1475.
1428
1429
<div id="cite-2"></div>
1430
'''[[#citeF-2|[2]]]'''   Aubry R, Idelsohn SR, Oñate E (2006)  Fractional step like schemes for free surface problems with thermal coupling using the Lagrangian PFEM. Computational Mechanics,  Vol. 38 (4-5), pp. 294&#8211;309.
1431
1432
<div id="cite-3"></div>
1433
'''[[#citeF-3|[3]]]'''  Belytschko T, Liu WK, Moran B, Non linear finite element for continua and structures, 2d Edition, Wiley, 2013.
1434
1435
<div id="cite-4"></div>
1436
'''[[#citeF-4|[4]]]'''  Carbonell JM, Oñate E, Suárez B (2010) Modeling of ground excavation with the Particle Finite Element Method. Journal of Engineering Mechanics (ASCE) 136(4):455&#8211;463
1437
1438
<div id="cite-5"></div>
1439
'''[[#citeF-5|[5]]]'''  Carbonell JM, Oñate E, Suárez B (2013) Modelling of tunnelling processes and cutting tool wear with the Particle Finite Element Method (PFEM). Comput. Mech. 52(3):607&#8211;629
1440
1441
<div id="cite-6"></div>
1442
'''[[#citeF-6|[6]]]'''  Chenot JL, Oñate E (1988) ''Modelling of Metal Forming Processes'', Kluwer Academic, Dordrecht.
1443
1444
<div id="cite-7"></div>
1445
'''[[#citeF-7|[7]]]'''   Donea J,  Huerta A  (2003) Finite element method for flow problems. J. Wiley
1446
1447
<div id="cite-8"></div>
1448
'''[[#citeF-8|[8]]]'''  Edelsbrunner H,  Mucke EP (1999) Three dimensional alpha shapes. ACM   Trans. Graphics  13:43&#8211;72
1449
1450
<div id="cite-9"></div>
1451
'''[[#citeF-9|[9]]]'''  Franci A, Oñate E, Carbonell JM (2013) On the effect of the  tangent bulk stiffness matrix in the analysis of free surface Lagrangian flows using PFEM. Research Report CIMNE PI402. Submitted to Int. J. Numer. Meth. Biomed. Engng.
1452
1453
<div id="cite-10"></div>
1454
'''[[#citeF-10|[10]]]'''  Franci A, Oñate E, Carbonell JM (2014) Unified updated Lagrangian formulation for fluid-structure interaction problems. Publication CIMNE No. PI404.
1455
1456
<div id="cite-11"></div>
1457
'''[[#citeF-11|[11]]]'''  Ghosh S, Castro JC,   Lee JK (2004) Material Processing and Design: Simulation and Applications. NUMIFORM 2004, American Institute of Physics.
1458
1459
<div id="cite-12"></div>
1460
'''[[#citeF-12|[12]]]'''  Huetnik J, Baaijens FPT  (Eds) (1998) ''Simulation of Materials Processing: Theory,  Methods and Applications, NUMIFORM 98''.  Enschede, Países Bajos,  A.A. Balkema, Rotterdam, 22&#8211;5 Jun.
1461
1462
<div id="cite-13"></div>
1463
'''[[#citeF-13|[13]]]'''  Idelsohn SR, Calvo N, Oñate E (2003c) Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng. 192(22-24):2649&#8211;2668
1464
1465
<div id="cite-14"></div>
1466
'''[[#citeF-14|[14]]]'''  Idelsohn SR, Oñate E, Del Pin F (2004)  The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int. J. Numer. Meth. Biomed. Engng.,  61(7):964&#8211;989
1467
1468
<div id="cite-15"></div>
1469
'''[[#citeF-15|[15]]]'''  Idelsohn SR, Marti J, Limache A, Oñate E (2008) Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluid-structure interaction problems via the PFEM. Comput Methods Appl Mech Engrg.  197:1762&#8211;1776
1470
1471
<div id="cite-16"></div>
1472
'''[[#citeF-16|[16]]]'''  Idelsohn SR, Mier-Torrecilla M, Oñate E  (2009) Multi-fluid flows with the Particle Finite Element Method. Comput Methods Appl Mech Engrg. 198:2750&#8211;2767
1473
1474
<div id="cite-17"></div>
1475
'''[[#citeF-17|[17]]]'''  Larese A,  Rossi R, Oñate E, Idelsohn SR (2008)  Validation of the Particle Finite Element Method (PFEM) for simulation of free surface flows. Engineering Computations 25(4):385&#8211;425
1476
1477
<div id="cite-18"></div>
1478
'''[[#citeF-18|[18]]]'''  Limache A, Idelsohn, SR, Rossi R, Oñate E (2007) The violation of objectivity in Laplace formulation of the Navier-Stokes equations. Int. J. Numerical Methods in Fluids, 54:639&#8211;664.
1479
1480
<div id="cite-19"></div>
1481
'''[[#citeF-19|[19]]]'''  Mori K (2001) ''Simulation of Material Processing: Theory, Methods and Applications'', A.A. Balkema, Lisse.
1482
1483
<div id="cite-20"></div>
1484
'''[[#citeF-20|[20]]]'''  Oliver X, Cante JC, Weyler R, González C, Hernández J (2007) Particle finite element methods in solid mechanics problems. In: Oñate E, Owen R (Eds) Computational Plasticity. Springer, Berlin, pp 87&#8211;103
1485
1486
<div id="cite-21"></div>
1487
'''[[#citeF-21|[21]]]'''  Oñate E (1998) Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Engng. 151:233&#8211;267
1488
1489
<div id="cite-22"></div>
1490
'''[[#citeF-22|[22]]]'''  Oñate E, Manzan M (1999)  A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. International Journal for Numerical Methods in Fluids,  Vol. 31, pp. 203&#8211;221.
1491
1492
<div id="cite-23"></div>
1493
'''[[#citeF-23|[23]]]'''  Oñate E (2000) A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comput Methods Appl Mech Engrg. 182(1&#8211;2):355&#8211;370
1494
1495
<div id="cite-24"></div>
1496
'''[[#citeF-24|[24]]]'''   Oñate E (2004) Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng. 60(1):255&#8211;281
1497
1498
<div id="cite-25"></div>
1499
'''[[#citeF-25|[25]]]'''   Oñate E, Rojek J, Taylor R, Zienkiewicz O (2004a) Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. Int. J. Numer. Meth. Engng. 59(11):1473&#8211;1500
1500
1501
<div id="cite-26"></div>
1502
'''[[#citeF-26|[26]]]'''  Oñate E, Idelsohn SR, Del Pin F, Aubry R (2004b) The particle   finite element method. An overview. Int. J. Comput. Methods    1(2):267&#8211;307
1503
1504
<div id="cite-27"></div>
1505
'''[[#citeF-27|[27]]]'''  Oñate E, M.A. Celigueta, Idelsohn SR (2006) Modeling bed erosion in   free surface flows by the Particle Finite Element Method, Acta   Geotechnia 1(4):237&#8211;252
1506
1507
<div id="cite-28"></div>
1508
'''[[#citeF-28|[28]]]'''  Oñate E, Valls A, García J  (2006) FIC/FEM formulation   with matrix stabilizing terms for incompressible flows at low and high   Reynold's numbers. Computational Mechanics 38 (4-5):440&#8211;455
1509
1510
<div id="cite-29"></div>
1511
'''[[#citeF-29|[29]]]'''  Oñate E, García J,  SR Idelsohn, Del Pin F (2006)  FIC   formulations for finite element analysis of incompressible flows. Eulerian,   ALE and Lagrangian approaches.  Comput. Meth. Appl. Mech. Engng.   195(23-24):3001&#8211;3037
1512
1513
<div id="cite-30"></div>
1514
'''[[#citeF-30|[30]]]'''  Oñate E, Rojek J, Chiumenti M, Idelsohn SR, Del Pin F, Aubry R (2006) Advances in stabilized finite element and particle methods for bulk forming processes. Comput. Meth. Appl. Mech. Engng. 195:6750&#8211;6777
1515
1516
<div id="cite-31"></div>
1517
'''[[#citeF-31|[31]]]'''  Oñate E, Owen DRJ, (Eds.), ''Computational Plasticity'', Springer 2007
1518
1519
<div id="cite-32"></div>
1520
'''[[#citeF-32|[32]]]'''  Oñate E, Idelsohn SR,  Celigueta MA, Rossi R (2008) Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. Comput. Meth. Appl. Mech. Engng. 197(19-20):1777&#8211;1800
1521
1522
<div id="cite-33"></div>
1523
'''[[#citeF-33|[33]]]'''  Oñate E (2009)  Structural analysis with the finite element method. Linear statics. Volume 1. Basis and Solids.  CIMNE-Springer
1524
1525
<div id="cite-34"></div>
1526
'''[[#citeF-34|[34]]]'''  Oñate E, Rossi R,  Idelsohn SR, Butler K (2010) Melting and spread of polymers in fire with the particle finite element method. Int. J. Numerical Methods in Engineering, 81(8):1046&#8211;1072
1527
1528
<div id="cite-35"></div>
1529
'''[[#citeF-35|[35]]]'''  Oñate E, Celigueta MA, Idelsohn SR,  Salazar F, Suárez B (2011) Possibilities of the particle finite element method for fluid-soil-structure interaction problems. Computational Mechanics, 48(3):307&#8211;318
1530
1531
<div id="cite-36"></div>
1532
'''[[#citeF-36|[36]]]'''  Oñate E, Idelsohn SR, Felippa C (2011) Consistent pressure Laplacian stabilization for incompressible continua via higher-order finite calculus. Int. J. Numer. Meth. Engng, 87 (1-5): 171&#8211;195
1533
1534
<div id="cite-37"></div>
1535
'''[[#citeF-37|[37]]]'''  Oñate E, Carbonell JM  (2013) Updated Lagrangian finite element formulation for quasi-incompressible fluids. Research Report PI393, CIMNE. Submitted to Comput. Mechanics
1536
1537
<div id="cite-38"></div>
1538
'''[[#citeF-38|[38]]]'''  Oñate E, Franci A, Carbonell JM (2013) Lagrangian formulation  for finite element analysis of quasi-incompressible fluids with reduced mass losses. Int. J. Num. Meth. Fluids, DOI: 10.1002/fld.3870
1539
1540
<div id="cite-39"></div>
1541
'''[[#citeF-39|[39]]]'''  Pietrzyk M, Kusiak J, Majta J, Hartley P, Pillinger I (2000) ''Metal Forming 2000''. A.A. Balkema, Rotterdam
1542
1543
<div id="cite-40"></div>
1544
'''[[#citeF-40|[40]]]'''  Pittman JFT, Zienkiewicz OC,  Wood RD, Alexander JM (Eds) (1984) '' Numerical Analysis of Forming Processes''. Wiley, Chichester
1545
1546
<div id="cite-41"></div>
1547
'''[[#citeF-41|[41]]]'''  Ryzhakov P, Oñate E, Rossi R, Idelsohn SR (2012) Improving mass conservation in simulation of incompressible flows. Int. J. Numer. Meth. Engng.  90(12):1435&#8211;1451
1548
1549
<div id="cite-42"></div>
1550
'''[[#citeF-42|[42]]]'''    Simo JC, Miehe C (1992)   Associative coupled thermoplasticity at finite strains Formulation, numerical analysis and implementation.   ''Computer Methods in Applied Mechanics and Engineering'', 98:41-104
1551
1552
<div id="cite-43"></div>
1553
'''[[#citeF-43|[43]]]'''  Tang B, Li JF, Wang TS (2009) Some improvements on free surface simulation by the particle finite element method. ''Int. J. Num. Methods in Fluids'', 60 (9):1032–-1054
1554
1555
<div id="cite-44"></div>
1556
'''[[#citeF-44|[44]]]'''  Zienkiewicz OC, Jain PC, Oñate E (1978) Flow of solids during forming and extrusion: some aspects of numerical solutions. ''Int. J. Solids Struct.'', 14:15&#8211;38
1557
1558
<div id="cite-45"></div>
1559
'''[[#citeF-45|[45]]]'''  Zienkiewicz OC, Oñate E,  Heinrich JC (1981)  A general formulation for coupled thermal flow of metals using finite elements, ''Int. J. Num. Meth. Eng.'', 17:1497&#8211;514, 1981
1560
1561
<div id="cite-46"></div>
1562
'''[[#citeF-46|[46]]]'''  Zienkiewicz OC,  Taylor RL, Zhu JZ (2005) The finite element method. The basis.  Vol. 1, 6th Ed., Elsevier
1563
1564
<div id="cite-47"></div>
1565
'''[[#citeF-47|[47]]]'''  Zienkiewicz OC,  Taylor RL (2005) The finite element method for   solid and structural mechanics. Vol. 2,  6th Ed., Elsevier
1566
1567
<div id="cite-48"></div>
1568
'''[[#citeF-48|[48]]]'''  Zienkiewicz OC,  Taylor RL, Nithiarasu P (2005)  The finite element   method for fluid dynamics.  Vol. 3, 6th Ed.,  Elsevier
1569

Return to Onate et al 2014c.

Back to Top

Document information

Published on 01/01/2014

Licence: CC BY-NC-SA license

Document Score

0

Views 32
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?