You do not have permission to edit this page, for the following reason:

You are not allowed to execute the action you have requested.


You can view and copy the source of this page.

x
 
1
2
3
==Abstract==
4
5
A new computational procedure for analysis of the melting and flame spread of polymers under fire conditions is presented. The method, termed Particle Finite Element Method (PFEM), combines concepts from particle-based techniques with those of the standard finite element method (FEM). The key feature of the PFEM is the use of an updated Lagrangian description to model the motion of nodes (particles) in the thermoplastic material. Nodes are viewed as material points which can freely move and even separate from the main analysis domain representing, for instance, the effect of melting and dripping of polymer particles. A mesh connects the nodes defining the discretized domain where the governing equations are solved as in the standard FEM. An incremental iterative scheme for the solution of the nonlinear transient coupled thermal-flow problem, including loss of mass by gasification, is used. Examples of the possibilities of the PFEM for the modelling and simulation of the melting and flame spread of polymers under different fire conditions are described. Numerical results are compared with experimental data provided by NIST.
6
7
8
9
'''keywords''' ''Melting, dripping, polymers, Particle Finite Element Method, PFEM''
10
11
12
==1 INTRODUCTION==
13
14
Thermoplastic polymer objects, including mattresses, upholstered furniture, and molded objects such as electronic housings and automobile parts, respond to fire by melting and dripping onto the surface below.  The flow of polymer material affects heat and mass transport within the object, and the accumulating melt pool below the object extends the flaming zone and increases the overall rate of heat release <span id='citeF-1'></span>[[#cite-1|[1]]]&#8211;<span id='citeF-3'></span>[[#cite-3|[3]]].  If the fire from the object and the pool fire interact, the intensity of the fire is enhanced even further. The spread rate of the melt pool and its burning behavior (including whether it is even able to sustain ignition) are affected by the flooring material as well as by the properties of the melt.
15
16
Computer modeling and simulation of the melting flow and flame spread of thermoplastics is extremely complex involving fluid flow, heat transfer, material degradation, flame chemistry, surface tension, and complex material properties <span id='citeF-1'></span>[[#cite-1|[1]]].  In addition, the drastic changes in shape pose a severe challenge to traditional modelling methods.  Attempts to model melt flow of polymeric material in fire using the volume of fluid (VOF) method have encountered difficulties with numerical instabilities and excessive runtimes <span id='citeF-4'></span>[[#cite-4|[4]]].
17
18
The Particle Finite Element Method (PFEM) <span id='citeF-5'></span>[[#cite-5|[5]]]&#8211;<span id='citeF-11'></span>[[#cite-11|[11]]] presented in this work is a powerful technique for modelling and analysis of complex multidisciplinary problems in fluid and solid mechanics involving coupled thermal effects, fragmentation and separation of fluid particles and fluid-structure interaction effects, among others. The PFEM combines Lagrangian particle-based techniques with the advantages of the integral formulation of the finite element method [[#cite-15|[15]]]. A key advantage of the Lagrangian formulation in the PFEM is the elimination of the convective terms in the fluid flow and thermal equations, which favours the simplicity (and symmetry) of the formulation, as well as computational efficiency.
19
20
In the PFEM, the particles represent the nodes of a finite element mesh.  The particles can move freely according to the velocity field, transporting their momentum and physical properties.  A robust and efficient remeshing algorithm connects the nodes into a finite element grid for solution of the state variables in the new configuration.  The PFEM has been used to solve a variety of free surface, fluid-structure interaction, and multiphase problems, including  breaking waves in harbours, ship hydrodynamics, dam bursting and metal casting, among many others <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[5,6]]], <span id='citeF-8'></span>[[#cite-8|[8]]]&#8211;<span id='citeF-11'></span>[[#cite-11|[11]]].
21
22
This paper describes the key aspects of the PFEM for analysis of the melting, flow and flame spread of polymer objects including loss of mass by gasification. The underlying ideas of the paper follow the original work of the authors in this field as presented in [Buetal07,Roetal07]. The recent developments of the PFEM described in this paper demonstrate its potential for solving complex melting flow problems for two and three dimensional heated polymer objects accounting for frictional contact and self-contact situations.
23
24
The paper is structured as follows. In the next section the basis of the PFEM is summarized. The essential governing equations and an overview of the discretization procedure and the general solution scheme are given. The potential of the PFEM for the simulation of the melt flow and spread of thermoplastic objects in fire is shown in examples of its application for different heated samples,including a rectangular slab, a triangular slab, polymer melt spreading over a surface, a chair, and a candle. Numerical results for the rectangular slab and the spreading melt are compared with experimental data obtained at NIST <span id='citeF-4'></span>[[#cite-4|[4]]].
25
26
==2 THE PARTICLE FINITE ELEMENT METHOD. AN OVERVIEW==
27
28
===2.1 Basic steps of the PFEM===
29
30
In the PFEM the analysis domain is modeled with an ''updated Lagrangian formulation'' <span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-11'></span>[[#cite-6|[6,8,11]]]. The analysis domain can include solid and fluid subdomains. As an example, we can model one or several thermoplastic objects and the surrounding air as forming part of the analysis domain. All variables in the fluid and solid domains are assumed to be known in the current configuration at time <math display="inline">t</math>. The new set of variables in both domains is sought  in the next configuration at time <math display="inline">t+\triangle t</math>. The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. To do this, the nodes discretizing the analysis domain are treated as material particles whose motion is tracked during the transient solution. This is useful to model the separation of particles from a solid domain, such as in the dripping of melt particles from a thermoplastic object, and to follow their subsequent motion as individual particles with a known density, an initial acceleration and velocity, and subject to gravity forces. Every node is a material point and hence is characterized by the density of the polymer material. The mass of a given domain is obtained by integrating the density at the different material points over the domain.
31
32
In this paper, ''we will consider the analysis domain to be formed of polymer material'' only, i.e. the surrounding hot air is not included as part of the analysis domain. Melted drips of the polymer material separating from the main body of the heated object will be treated as independent analysis domains.
33
34
The way the PFEM solution process operates for the problems we are solving in this paper is schematically shown in Figure [[#img-1|1]]. The figure represents a polymer object hanging from a wall and subjected to an incoming heat flux <math display="inline">q</math> acting at the lower part of the object.
35
36
The ''collection or cloud of nodes'' pertaining to the polymer analysis domain will be defined as (<math display="inline">C</math>), the volume defining the thermoplastic analysis domain as (<math display="inline">V</math>), and the mesh discretizing this domain as (<math display="inline">M</math>).
37
38
A typical solution with the PFEM involves the following steps:
39
40
<ol>
41
42
<li>The starting point at each time step is the cloud of points in the polymer     domain. For instance, <math display="inline">^\circ C</math> and <math display="inline">^{n}C</math> are the clouds at the initial time and at time <math display="inline">t = t_n</math>, respectively (Figure [[#img-1|1]]). </li>
43
<li>Identify the boundaries defining the analysis domain <math display="inline">^{n}V</math>. This is an     essential step as some boundaries, such as the free surface in the melting     zone, may be severely distorted during the solution process. Some nodes may     even separate from the boundary.     The Alpha-Shape method <span id='citeF-14'></span>[[#cite-14|[14]]] is used for the boundary definition. Note that the     analysis domain can include different subdomains formed by the main polymer     object and the dripping and spread zones (Figure [[#img-1|1]]). </li>
44
<li>Discretize the polymer analysis domain with a finite element mesh <math display="inline">^{n}M</math>.     In our work an innovative mesh generation scheme based on the extended Delaunay     tessellation is typically used [[#cite-7|[7]]]. </li>
45
<li>Solve the coupled Lagrangian equations of motion for the polymer domains.     Compute the relevant state variables at the next (updated) configuration for <math display="inline">t+\triangle t</math>:     velocities, strain rates, strains, pressure, viscous stresses and temperature     in the polymer.  </li>
46
<li>Move the mesh nodes to a new position <math display="inline">^{n+1}C</math>, where <math display="inline">n+1</math> denotes the time <math display="inline">t_n+\triangle t</math>, in terms of the time increment size.  For the fully Lagrangian     application discussed here, the nodes are moved according to the forces determined in step 4.  </li>
47
<li>Go back to step 1 and repeat the solution process for the next time step. </li>
48
49
</ol>
50
51
<div id='img-1'></div>
52
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
53
|-
54
|[[Image:Draft_Samper_846581836-Figure1.png|600px|Polymer object subjected to a heat flux q applied to its lower boundary (arrows indicate incoming heat flux). Sequence of steps to update the “cloud” of nodes representing the object in time using the PFEM]]
55
|- style="text-align: center; font-size: 75%;"
56
| colspan="1" | '''Figure 1:''' Polymer object subjected to a heat flux <math>q</math> applied to its lower boundary (arrows indicate incoming heat flux). Sequence of steps to update the “cloud” of nodes representing the object in time using the PFEM
57
|}
58
59
It is important to note that the quality of the numerical solution depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can therefore be used to refine the mesh at the onset of the solution for each the step in order to improve the solution in zones where large motions of the fluid or the structure occur. Indeed, remeshing implies the introduction of new nodes in the mesh for which the state variables must be assigned via a projection method. The error introduced by the projection is corrected using an implicit algorithm which sought for a converged solution at the next configuration at time. In the examples shown in this paper the number of nodes during the transient solution has been kept constant.
60
61
Two distinct features of the PFEM are its capability to model the separation of material fragments from the main body of the object and its ability to model ''frictional contact situations''. The contact between two solid interfaces is treated by introducing a layer of ''contact elements'' between the two interacting domains. This layer is ''automatically created during the mesh generation step'' by prescribing a minimum distance (<math display="inline">h_{c}</math>) between the two interacting boundaries that controls the accuracy of the contact model. If the distance is greater than the minimum value (<math display="inline">h_{c}</math>) then the generated elements are given the properties of air and treated as standard fluid elements, or else they are removed from the analysis domain (Figure [[#img-2|2]]). Otherwise, the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacements (or velocities) is introduced so as to model frictional contact effects. This technique also prevents nodes in a fluid domain from crossing over rigid boundaries <span id='citeF-8'></span><span id='citeF-11'></span>[[#cite-8|[8,11]]].
62
63
The PFEM has proven to be very effective for modeling complex frictional contact conditions between two or more interacting rigid or deformable bodies or between fluid and solid interfaces in an extremely simple manner. Examples of applications of the PFEM can be found in <span id='citeF-6'></span>[[#cite-6|[6]]]&#8211;<span id='citeF-13'></span>[[#cite-13|[13]]].
64
65
<div id='img-2'></div>
66
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
67
|-
68
|[[Image:Draft_Samper_846581836-Fig2a.png|540px|]]
69
|[[Image:Draft_Samper_846581836-Fig2b.png|540px|Modelling of contact between the melting object and a fixed boundary]]
70
|- style="text-align: center; font-size: 75%;"
71
| colspan="2" | '''Figure 2:''' Modelling of contact between the melting object and a fixed boundary
72
|}
73
74
===2.2 Governing equations===
75
76
It is assumed that the polymer melt flow is governed by the equations of an incompressible fluid with a temperature dependent viscosity. A quasi-rigid behaviour of the polymer object at room temperature is reproduced by setting the viscosity to a sufficiently high value that the unheated polymer moves a negligible distance over the duration of the problem. As temperature increases in the thermoplastic object due to heat exposure, the viscosity decreases by several orders of magnitude as a function of temperature. This induces the melt and flow of the particles in the heated zone. The key equations to be solved in the polymer melt flow problem, ''written in the Lagrangian frame of reference'', are the following:
77
78
==Momentum==
79
80
<span id="eq-1"></span> 
81
82
{| class="formulaSCP" style="width: 100%; text-align: left;" 
83
|-
84
| 
85
{| style="text-align: left; margin:auto;width: 100%;" 
86
|-
87
| style="text-align: center;" | <math>\rho \frac{dv_i}{dt}= \frac{\partial \sigma _{ij}}{\partial x_j}+b_i \quad \hbox{in } \Omega  </math>
88
|}
89
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
90
|}
91
92
==Mass balance==
93
94
<span id="eq-2"></span> 
95
96
{| class="formulaSCP" style="width: 100%; text-align: left;" 
97
|-
98
| 
99
{| style="text-align: left; margin:auto;width: 100%;" 
100
|-
101
| style="text-align: center;" | <math>\frac{\partial v_{i}}{\partial x_i}=0 \quad \hbox{in } \Omega  </math>
102
|}
103
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
104
|}
105
106
==Heat transport==
107
108
<span id="eq-3"></span> 
109
110
{| class="formulaSCP" style="width: 100%; text-align: left;" 
111
|-
112
| 
113
{| style="text-align: left; margin:auto;width: 100%;" 
114
|-
115
| style="text-align: center;" | <math>\rho c \frac{dT}{dt}= \frac{\partial }{\partial x_i} \left(k_i  \frac{\partial T}{\partial x_i} \right) +Q \quad \hbox{in } \Omega  </math>
116
|}
117
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
118
|}
119
120
In the above equations <math display="inline">v_{i}</math> is the velocity along the ''i''th global (cartesian) axis, <math display="inline">T</math> is the temperature, <math display="inline">\rho </math>, <math display="inline">c</math> and <math display="inline">k_i</math> are the density (assumed constant), the specific heat and the conductivity of the material along the ''i''th coordinate direction, respectively, <math display="inline">b_i</math> and <math display="inline">Q</math> are the body forces and the heat source per unit volume, respectively, and <math display="inline">\sigma _{ij}</math> are the (Cauchy) stresses related to the velocities by the standard constitutive equation (for an incompressible Newtonian fluid)
121
122
<span id="eq-4"></span> 
123
124
{| class="formulaSCP" style="width: 100%; text-align: left;" 
125
|-
126
| 
127
{| style="text-align: left; margin:auto;width: 100%;" 
128
|-
129
| style="text-align: center;" | <math>\sigma _{ij} =s_{ij}-p\delta _{ij} </math>
130
|}
131
| style="width: 5px;text-align: right;white-space: nowrap;" | (4.a)
132
|}
133
134
{| class="formulaSCP" style="width: 100%; text-align: left;" 
135
|-
136
| 
137
{| style="text-align: left; margin:auto;width: 100%;" 
138
|-
139
| style="text-align: center;" | <math>s_{ij}=2\mu \varepsilon _{ij}\quad ,\quad \varepsilon _{ij}= \frac{1}{2}\left( {\partial v_i \over \partial x_j}+{\partial v_j \over \partial x_i} \right) </math>
140
|}
141
| style="width: 5px;text-align: right;white-space: nowrap;" | (4.b)
142
|}
143
144
In Eqs.([[#eq-4|4]]) <math display="inline">s_{ij}</math>  are the deviatoric stresses, <math display="inline">p</math>  is the pressure (assumed to be positive in compression), <math display="inline">\varepsilon _{ij}</math> is the rate of deformation and <math display="inline">\mu </math> is the viscosity. In the following we will assume <math display="inline">k_i=k</math> and the viscosity <math display="inline">\mu </math> to be a known function of temperature, i.e. <math display="inline">\mu =\mu (T)</math>.
145
146
In Eqs. ([[#eq-1|1]]) &#8211; ([[#eq-3|3]]) <math display="inline"> \frac{d\varphi }{dt}</math> means the total (material) derivative of any variable <math display="inline">\varphi </math>  with respect to time.
147
148
Indexes in Eqs.([[#eq-1|1]]) &#8211; ([[#eq-4|4]]) range from <math display="inline">i,j=1,n_d</math>, where <math display="inline">n_d</math> is the number of space dimensions of the problem (i.e. <math display="inline">n_d  = 2</math> for two-dimensional problems).
149
150
Eqs.([[#eq-1|1]]) &#8211; ([[#eq-4|4]]) are completed with the standard boundary conditions of prescribed velocities and surface tractions in the mechanical problem and prescribed temperature and prescribed normal heat flux in the thermal problem [[#cite-6|[6,8,9,10]]]. For surfaces exposed to fire conditions, energy losses due to radiation and convection must be taken into account, and the thermal boundary condition is:
151
152
<span id="eq-5.a"></span>
153
{| class="formulaSCP" style="width: 100%; text-align: left;" 
154
|-
155
| 
156
{| style="text-align: left; margin:auto;width: 100%;" 
157
|-
158
| style="text-align: center;" | <math>k\frac{\partial T}{\partial n}+\bar q_n =0 </math>
159
|}
160
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.a)
161
|}
162
163
with
164
165
<span id="eq-5.b"></span>
166
{| class="formulaSCP" style="width: 100%; text-align: left;" 
167
|-
168
| 
169
{| style="text-align: left; margin:auto;width: 100%;" 
170
|-
171
| style="text-align: center;" | <math>\bar q_n =q_n +\epsilon \sigma (T^4-T^4_0) + \alpha _c(T-T_{0}) \quad \hbox{in }\Gamma _q </math>
172
|}
173
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.b)
174
|}
175
176
where <math display="inline">\frac{\partial T}{\partial n}</math> is the derivative of temperature within the object along the direction normal to the Neumann boundary <math display="inline">\Gamma _q</math>, <math display="inline">q_n</math> is the outgoing heat flux per unit area along the normal to the boundary surface, and the last two terms of Eqs.([[#eq-5.b|5.b]]) are the radiative and convective heat losses respectively, where <math display="inline">T_0</math> is the reference temperature, <math display="inline">\epsilon </math> is the emissivity, <math display="inline">\sigma </math> is the Stefan-Boltzmann constant, and <math display="inline">\alpha _c</math> is the convection coefficient.
177
178
We note that Eqs.([[#eq-1|1]]) &#8211; ([[#eq-4|4]]) are the standard ones for modeling the deformation of viscoplastic materials using the so-called “flow approach” <span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-19|[19,20]]]. In our work we have limited the viscosity to be a non linear function of the temperature. Other dependencies of the viscosity with the shear rate and the yield stress, typical of the non-Newtonian flow of solids, can be easily incorporated in the model, thus broadening  its applicability to other materials.
179
180
===2.3 Discretization of the equations===
181
182
A key problem in the numerical solution of ([[#eq-1|1]]) &#8211; ([[#eq-4|4]]) is the satisfaction of the incompressibility condition (Eq.[[#eq-2|(2)]]). A number of procedures to solve this problem exist in the finite element literature <span id='citeF-11'></span>[[#cite-11|[11]]]. In our approach, we use a stabilized formulation based on the so-called finite calculus procedure <span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span>[[#cite-16|[16,17,18]]]. The essence of this method is the solution of a ''modified mass balance'' equation that is written as
183
184
<span id="eq-6"></span> 
185
186
{| class="formulaSCP" style="width: 100%; text-align: left;" 
187
|-
188
| 
189
{| style="text-align: left; margin:auto;width: 100%;" 
190
|-
191
| style="text-align: center;" | <math>{\partial v_i \over \partial x_i} + \sum \limits _{i=1}^3 \tau \left[{\partial p \over \partial x_i} +\pi _i\right]=0 </math>
192
|}
193
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
194
|}
195
196
where <math display="inline">\tau </math> is a stabilization parameter given by <span id='citeF-17'></span>[[#cite-17|[17]]]
197
198
<span id="eq-7"></span> 
199
200
{| class="formulaSCP" style="width: 100%; text-align: left;" 
201
|-
202
| 
203
{| style="text-align: left; margin:auto;width: 100%;" 
204
|-
205
| style="text-align: center;" | <math>\tau = \left(\frac{2\rho \vert \mathbf{v}\vert }{h}+\frac{8\mu }{3h^2} \right)^{-1} </math>
206
|}
207
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
208
|}
209
210
In the above, <math display="inline">h</math> is a characteristic length of each finite element (such as <math display="inline">[A^{(e)}]^{1/2}</math> for 2D elements), and <math display="inline">\vert \mathbf{v}\vert </math> is the modulus of the velocity vector. Also in Eq.[[#eq-5.a|(5)]] <math display="inline">\pi _i</math> are auxiliary pressure projection variables chosen so as to ensure that in the second term in Eq.[[#eq-6|(6)]] <math display="inline">\pi _i</math> can be interpreted as a weighted sum of the residuals of the momentum equations and therefore it vanishes for the exact solution. The set of governing equations for the velocities, the pressure and the <math display="inline">\pi _i</math> variables is completed by adding the following constraint equation to the set of governing equations
211
212
<span id="eq-8"></span> 
213
214
{| class="formulaSCP" style="width: 100%; text-align: left;" 
215
|-
216
| 
217
{| style="text-align: left; margin:auto;width: 100%;" 
218
|-
219
| style="text-align: center;" | <math>\int _V \tau w_i\left({\partial p \over \partial x_i}  +\pi _i\right)dV =0 \quad ,\quad i=1,n_d </math>
220
|}
221
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
222
|}
223
224
where <math display="inline">w_i</math> are arbitrary weighting functions.
225
226
The rest of the integral equations are obtained by applying the standard Galerkin technique to the governing equations [[#eq-1|(1)]], [[#eq-2|(2)]] and [[#eq-3|(3)]]  and the appropriate boundary conditions for the numerical and thermal problems<span id='citeF-8'></span><span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-8|[8,10,11]]].
227
228
We interpolate in the standard finite element fashion the set of problem variables. For 3D problems these are the three velocities <math display="inline">v_i</math>, the pressure <math display="inline">p</math>, the temperature <math display="inline">T</math> and the three pressure gradient projections <math display="inline">\pi _i</math>. In our work we use an equal order ''linear interpolation'' for all variables over meshes of 3-noded triangles (in 2D) and 4-noded tetrahedra (in 3D) <span id='citeF-8'></span><span id='citeF-11'></span><span id='citeF-15'></span>[[#cite-8|[8,11,15]]]. The resulting set of discretized equations has the following form
229
230
==Momentum==
231
232
<span id="eq-9"></span> 
233
234
{| class="formulaSCP" style="width: 100%; text-align: left;" 
235
|-
236
| 
237
{| style="text-align: left; margin:auto;width: 100%;" 
238
|-
239
| style="text-align: center;" | <math>{\boldsymbol M}\dot{\bar{\boldsymbol v}} + {\boldsymbol K}(\mu ) \bar {\boldsymbol v} - {\boldsymbol G} \bar {\boldsymbol p}= {\boldsymbol  f} </math>
240
|}
241
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
242
|}
243
244
==Mass balance==
245
246
<span id="eq-10"></span> 
247
248
{| class="formulaSCP" style="width: 100%; text-align: left;" 
249
|-
250
| 
251
{| style="text-align: left; margin:auto;width: 100%;" 
252
|-
253
| style="text-align: center;" | <math>{\boldsymbol G}^T \bar {\boldsymbol v} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}^T\bar {\boldsymbol    \pi }={\boldsymbol 0} </math>
254
|}
255
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
256
|}
257
258
==Pressure gradient projections==
259
260
<span id="eq-11"></span> 
261
262
{| class="formulaSCP" style="width: 100%; text-align: left;" 
263
|-
264
| 
265
{| style="text-align: left; margin:auto;width: 100%;" 
266
|-
267
| style="text-align: center;" | <math>\hat {\boldsymbol M} \bar {\boldsymbol \pi } + {\boldsymbol Q}^T \bar {\boldsymbol p} ={\boldsymbol 0} </math>
268
|}
269
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
270
|}
271
272
==Heat transfer==
273
274
<span id="eq-12"></span> 
275
276
{| class="formulaSCP" style="width: 100%; text-align: left;" 
277
|-
278
| 
279
{| style="text-align: left; margin:auto;width: 100%;" 
280
|-
281
| style="text-align: center;" | <math>{\boldsymbol C}\dot{\bar{\boldsymbol T}} +{\boldsymbol H}\bar {\boldsymbol T}= {\boldsymbol q} </math>
282
|}
283
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
284
|}
285
286
In Eqs.([[#eq-9|9]]) &#8211; ([[#eq-12|12]]) <math display="inline">\bar {(\cdot )}</math> denotes nodal variables, <math display="inline">\dot{\bar {(\cdot  )}}= \frac{d}{dt}  \bar {(\cdot )}</math>. The different matrices and vectors are given in the Appendix.
287
288
The solution in time of Eqs.([[#eq-9|9]]) &#8211; ([[#eq-12|12]]) can be performed using any time integration scheme typical of the updated Lagrangian finite element method. A basic algorithm following the conceptual process described in Section [[#2.1 Basic steps of the PFEM|2.1]] is presented in Box I <span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-15'></span>[[#cite-6|[6,8,9,10,15]]]. In Box I <math display="inline">^{t+1} \bar{\boldsymbol a}^{i+1}</math>  denotes the values of the nodal variables <math display="inline">\bar{\boldsymbol a}</math> at time <math display="inline">t+\Delta t</math> and iteration <math display="inline">i+1</math>. Note that the position of the analysis domain in the next equilibrium configuration at time <math display="inline">t+\Delta  t</math> is also an output of the solution, as is typical in updated Lagrangian methods. We also recall the coupling of the flow and thermal equations via the dependence of the viscosity <math display="inline">\mu </math> with the temperature. The convergence of the solution process at each time step ensures the stability and consistency of the formulation.
289
Within a time step (once the discretization is fixed) the mass and volume are conserved in a variational sense.  Upon remeshing, the volume is not conserved in a strict sense, but the alpha-shape method  provides such conservation in a statistical sense <span id='citeF-6'></span><span id='citeF-8'></span>[[#cite-6|[6,8]]].
290
291
==3 ACCOUNTING FOR GASIFICATION EFFECTS==
292
293
The effect of gasification can be introduced by adding a (nonlinear) energy loss term in the energy equation. This term represents the energy that migrates from the condensed phase object to the gas due to the gasification of a part of the material during the heating process. The volumetric gasification heat flux has the following form
294
295
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
296
|-
297
298
1. LOOP OVER TIME STEPS, <math>t=1</math>, NTIME
299
300
Known values
301
302
<math>{~}^t\bar {\boldsymbol x},{~}^t\bar {\boldsymbol v},{~}^t\bar {\boldsymbol p},{~}^t\bar{\boldsymbol \pi }, {~}^t\bar {\boldsymbol T},{~}^t\mu , {~}^t{\boldsymbol f}, {~}^t{\boldsymbol q}, {~}^tC, {~}^tV, {~}^tM</math>
303
304
305
2. LOOP OVER NUMBER OF ITERATIONS, <math>i=1</math>, NITER
306
307
<math>\bullet </math> Compute the nodal velocities by solving Eq. [[#eq-9|(9)]]
308
309
<math>\displaystyle \left[\frac{1}{\Delta t}{\boldsymbol M} +{\boldsymbol K}\right]{~}^{t+1}\bar  {\boldsymbol v}^{i+1}= {~}^{t+1}{\boldsymbol f}+{\boldsymbol G} {~}^{t+1}\bar {\boldsymbol p}^i +\frac{1}{\Delta  t}{\boldsymbol M} {~}^t\bar {\boldsymbol v}</math>
310
311
<math>\bullet </math> Compute  nodal pressures from Eq.[[#eq-10|(10)]]
312
313
<math>{\boldsymbol L} {~}^{t+1}\bar {\boldsymbol p}^{i+1}=-{\boldsymbol G}^T {~}^{t+1}\bar {\boldsymbol v}^{i+1}- {\boldsymbol  Q}{~}^{t+1}\bar{\boldsymbol \pi }^i</math>
314
315
<math>\bullet </math> Compute nodal pressure gradient projections from Eq. [[#eq-11|(11)]]
316
317
<math>{~}^{t+1}\bar {\boldsymbol \pi }^{i+1}= - \hat {\boldsymbol M}_D [{\boldsymbol Q}^T] {~}^{t+1}\bar {\boldsymbol  p}^{i+1}\quad ,\quad \hat {\boldsymbol M}_D=diag[\hat {\boldsymbol M}]</math>
318
319
<math>\bullet </math> Compute nodal temperatures from Eq.[[#eq-12|(12)]]
320
321
<math>\displaystyle \left[\frac{1}{\Delta t}{\boldsymbol C} +{\boldsymbol H}\right]{~}^{t+1}\bar  {\boldsymbol T}^{i+1}={~}^{t+1}{\boldsymbol q}-\frac{1}{\Delta t}{\boldsymbol C} {~}^t\bar {\boldsymbol T}</math>
322
323
<math>\bullet </math> Update position of analysis domain nodes:
324
325
<math>{~}^{t+1}\bar {\boldsymbol x}^{i+1}= {~}^t{\boldsymbol x}^i+ {~}^{t+1}{\boldsymbol  v}^{i+1}\Delta t</math>
326
327
328
Define new “cloud” of nodes <math>^{t+1} {C}^{i+1}</math>
329
330
<math>\bullet </math> Update viscosity values in terms of temperature
331
332
<math>{~}^{t+1}\mu =\mu (^{t+1}\bar {\boldsymbol T}^{i+1})</math>
333
334
335
Check convergence  for all variables <math>\to </math> NO <math>\to </math>  Next iteration <math>i\to i+1</math>
336
337
<math>~~~~~~~~~~\downarrow </math> YES
338
339
340
Next time step <math>t\to t+1</math>
341
342
<math>\bullet </math> Identify new analysis domain boundary: <math>^{t+1}V</math>
343
344
<math>\bullet </math> Generate mesh: <math>{~}^{t+1}M</math>
345
346
Go to 1
347
348
'''Box I'''. Flow chart of basic PFEM algorithm for the polymer domain
349
350
|}
351
352
353
354
<span id="eq-13"></span> 
355
356
{| class="formulaSCP" style="width: 100%; text-align: left;" 
357
|-
358
| 
359
{| style="text-align: left; margin:auto;width: 100%;" 
360
|-
361
| style="text-align: center;" | <math>q_{gas}=\rho H\bar \varepsilon _v </math>
362
|}
363
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
364
|}
365
366
where <math display="inline">H</math> is the enthalpy of vaporization,
367
368
<span id="eq-14"></span> 
369
370
{| class="formulaSCP" style="width: 100%; text-align: left;" 
371
|-
372
| 
373
{| style="text-align: left; margin:auto;width: 100%;" 
374
|-
375
| style="text-align: center;" | <math>\bar \varepsilon _v =f(T) </math>
376
|}
377
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
378
|}
379
380
where <math display="inline">f(T)</math> expresses the (generally nonlinear)  relation between the rate of volume variation due to the temperature, <math display="inline">\bar   \varepsilon _v</math>,  and the temperature itself. In our work, the following Arrhenius function is  chosen, consistent with single-step, first-order thermal decomposition kinetics <span id='citeF-4'></span>[[#cite-4|[4]]]
381
382
<span id="eq-15"></span> 
383
384
{| class="formulaSCP" style="width: 100%; text-align: left;" 
385
|-
386
| 
387
{| style="text-align: left; margin:auto;width: 100%;" 
388
|-
389
| style="text-align: center;" | <math>f(T) = - Ae^{-\frac{E}{RT}}  </math>
390
|}
391
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
392
|}
393
394
where <math display="inline">A</math> is the pre-exponential function, <math display="inline">E</math> is the activation energy, <math display="inline">R</math> is the universal gas constant, and <math display="inline">T</math> is the absolute temperature expressed in Kelvin.
395
396
Once the temperature field is known at each iteration of a time step, the rate of volume variation <math display="inline">\bar \varepsilon _v</math> is fixed at every point of the mesh. This allows defining a continuum distribution of the term <math display="inline">\bar \varepsilon _v</math>   over the whole polymer domain.
397
398
The computed mass loss rate has to be included in the problem to ensure that the volume variation of the sample is correctly modeled. The approach is thus to solve the momentum and heat transfer equations, prescribing as constraints the local rate of volume variation and the gasification heat flux.
399
400
From a practical point of view this simply implies adding the gasification heat flux as an additional volumetric heat source term in vector <math display="inline">\mathbf{q}</math> in the heat transfer equations (see [[#APPENDIX|Appendix]]) and adding the volumetric deformation rate <math display="inline">\bar \varepsilon _v</math> into the stabilized mass balance equation ([[#eq-6|6]]) as
401
402
<span id="eq-16"></span> 
403
404
{| class="formulaSCP" style="width: 100%; text-align: left;" 
405
|-
406
| 
407
{| style="text-align: left; margin:auto;width: 100%;" 
408
|-
409
| style="text-align: center;" | <math>{\partial v_i \over \partial x_i}-\bar \varepsilon _v +\sum \limits _{i=1}^3 \tau _i \left[{\partial p \over \partial x_i} +\pi _i\right]=0 </math>
410
|}
411
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
412
|}
413
414
The new discretized mass balance equation is
415
416
<span id="eq-17"></span> 
417
418
{| class="formulaSCP" style="width: 100%; text-align: left;" 
419
|-
420
| 
421
{| style="text-align: left; margin:auto;width: 100%;" 
422
|-
423
| style="text-align: center;" | <math>{\boldsymbol G}^T \bar  {\boldsymbol v} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q} \bar {\boldsymbol \pi }={\boldsymbol f}_{g} </math>
424
|}
425
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
426
|}
427
428
where <math display="inline">{\boldsymbol f}_g</math> with <math display="inline">f_{g_i}= \int _{V^e} N_i \bar \varepsilon _v dV</math> is the forcing term contributed by the gasification heat flux.
429
430
An example of the effects of gasification will be shown in the following section.
431
432
==4 NUMERICAL EXAMPLES==
433
434
The examples were run in the PFIRE code in which the formulation described above has been implemented using the Kratos software platform <span id='citeF-23'></span>[[#cite-23|[23]]].
435
436
===4.1 Validation of the gasification model===
437
438
A simple computational test  has been chosen to verify the accuracy of the gasification algorithm.
439
440
A <math display="inline">5 \times 10</math> cms horizontal slab is insulated along its sides and bottom by an impermeable wall and is subjected to a heat flux of 30 kW/m<math display="inline">^2</math> acting on the free boundary. The properties of the material are shown in Table I. The analysis was carried out with several meshes of 3-noded triangles.
441
442
The example is approximately equivalent to a 1D gasification model for which the possibility exist of finding an accurate semi-analytical solution using Mathematica. The horizontal setup guarantees that the only mass loss is induced by the gasification term as no flow is allowed to leave the sample due to the presence of the impermeable boundaries <span id='citeF-12'></span>[[#cite-12|[12]]].
443
444
The results of  mass versus time for the different meshes considered are represented in Figure [[#img-3|3]].
445
446
<div id='img-3'></div>
447
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
448
|-
449
|[[Image:Draft_Samper_846581836-Fig3_con293.png|600px|Gasification of a rectangular slab (5×10 cms)   Analytical (1D model) and PFEM results for different mesh sizes.]]
450
|- style="text-align: center; font-size: 75%;"
451
| colspan="1" | '''Figure 3:''' Gasification of a rectangular slab (<math>5\times 10</math> cms)   Analytical (1D model) and PFEM results for different mesh sizes.
452
|}
453
454
Results for the finer mesh show a sufficiently close agreement between the numerical results and the analytical 1D solution, confirming the correctness of the gasification approach used.
455
456
Some screenshots of the sample geometry at different times  for the finer mesh are shown in Figure [[#img-4|4]].
457
458
Although the results are satisfactory, it is important to observe that a high mesh density is required in the vicinity of the free surface to capture accurately the analytical solution. This is due to the sharp gradient exhibited by the Arrhenius function in the proximity of the free surface, which needs to be accurately resolved. Such situation calls for an adaptive mesh refinement algorithm. This possibility will be implemented in forthcoming works.
459
460
<div id='img-4'></div>
461
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
462
|-
463
|[[Image:Draft_Samper_846581836-malla_gasification.png|360px|]]
464
|[[Image:Draft_Samper_846581836-temp_gasification_2000.png|420px|]]
465
|-
466
|[[Image:Draft_Samper_846581836-temp_gasification_4000.png|420px|]]
467
|[[Image:Draft_Samper_846581836-temp_gasification_6000.png|420px|Gasification of a heated slab. Initial geometry and mesh and slab shapes at 2000s, 4000s   and 6000s after onset of gasification.]]
468
|- style="text-align: center; font-size: 75%;"
469
| colspan="2" | '''Figure 4:''' Gasification of a heated slab. Initial geometry and mesh and slab shapes at 2000s, 4000s   and 6000s after onset of gasification.
470
|}
471
472
===4.2 Melting and flow of a rectangular slab===
473
474
In the next example shown, the PFEM is used to simulate an experiment performed at NIST in which a slab of polymeric material is mounted vertically and exposed to uniform radiant heating on one face.  Degradation of the polymer decreases its viscosity by several orders of magnitude and produces fuel gases. Polymer melt is captured by a pan below the sample.
475
476
A schematic of the apparatus used in the experiments is shown in Figure [[#img-5|5]]. A rectangular polymeric sample of dimensions 10 cm high by 10 cm wide by 5 cm thick is mounted upright and exposed to uniform heating on one face from a radiant cone heater placed on its side.   The sample is insulated on its lateral and rear faces.  The melt flows down the heated face of the sample and drips onto a surface below.  A load cell monitors the mass of polymer remaining in the sample, and a laboratory balance measures the mass of polymer falling onto the catch surface. Details of the experimental setup are given in previous publications [BuOhLi04,Buetal07,Roetal07].
477
478
<div id='img-5'></div>
479
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
480
|-
481
|[[Image:Draft_Samper_846581836-Fig3_polymers.png|338px|Polymer melt experiment. Viscosity vs. temperature for PP702N   polypropylene in its initial undegraded form and after exposure to 30   kW/m² and 40 kW/m² heat fluxes.  The black curve follows the   extrapolation of viscosity to high temperatures.]]
482
|- style="text-align: center; font-size: 75%;"
483
| colspan="1" | '''Figure 5:''' Polymer melt experiment. Viscosity vs. temperature for PP702N   polypropylene in its initial undegraded form and after exposure to 30   kW/m<math>^{2}</math> and 40 kW/m<math>^{2}</math> heat fluxes.  The black curve follows the   extrapolation of viscosity to high temperatures.
484
|}
485
486
487
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
488
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Parameter values for polypropylene type PP702N taken from <span id='citeF-22'></span>[[#cite-22|[22]]]. The value of <math>\epsilon =1.0</math> for the emissivity has been assumed
489
|-
490
| style="text-align: left;" |    Parameter 
491
| Symbol 
492
| Value 
493
| Units
494
|-
495
| style="text-align: left;" |    Density 
496
| <math>\rho </math>
497
| 900 
498
| kg/m<math display="inline">^3</math> 
499
|-
500
| style="text-align: left;" | Thermal conductivity 
501
| <math>k</math>
502
| 0.25 
503
| W/m-K  
504
|-
505
| style="text-align: left;" | Specific heat 
506
| <math>c</math>
507
| 2400 
508
| J/kg-K  
509
|-
510
| style="text-align: left;" | Reference temperature 
511
| <math>T_0</math>
512
| 298 
513
| K  
514
|-
515
| style="text-align: left;" | Emissivity 
516
| <math>\epsilon </math>
517
| 1.0 
518
| - 
519
|-
520
| style="text-align: left;" | Convective heat transfer coefficient 
521
| <math>\alpha _c</math>
522
| 8 
523
| W/m<math display="inline">^2</math>-K 
524
|-
525
| style="text-align: left;" | Heat of vaporization 
526
| <math>H</math>
527
| <math>8\times 10^5</math>
528
| J/kg  
529
|-
530
| style="text-align: left;" | Pre-exponential function 
531
| <math>A</math>
532
| <math>2.18\times 10^{12}</math>
533
| s<math display="inline">^{-1}</math> 
534
|-
535
| style="text-align: left;" | Activation energy divided by universal
536
|-
537
| style="text-align: left;" | gas constant 
538
| <math>E/R</math>
539
| 24400 
540
| K 
541
|-
542
| style="text-align: left;" | Stefan-Boltzmann constant 
543
| <math>\sigma </math>
544
| 5.67<math display="inline">\times 10^{-8}</math>
545
| W/m<math display="inline">^{2}</math>- K<math display="inline">^{4}</math>
546
547
|}
548
549
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
550
|-
551
552
<pre>
553
554
Tc: Temperature in degrees Celsius; T: Temperature in Kelvins
555
556
def f1(Tc):
557
return 10**(14.48 - 0.13858*Tc + 5.5960e-4*Tc*Tc - 7.8665e-7*Tc*Tc*Tc)
558
559
def f2(Tc):
560
return 10**(53.19 - 0.2542*Tc + 2.9879e-4*Tc*Tc)
561
562
def AuxFunction(T):
563
564
Tc = T - 273
565
if( Tc <= 25):
566
mu = 1e6
567
elif (Tc > 25 and  Tc <= 200):
568
mu = 1e6*(200-Tc)/175 + f1(200)
569
elif(Tc > 200 and Tc <= 350):
570
mu = f1(Tc)
571
elif(Tc > 350 and Tc < 425 ):
572
mu = f2(Tc)
573
else:
574
mu = f2(425)
575
576
return mu
577
</pre>
578
579
580
|}
581
582
The first computational representation of this problem is a 2-D rectangle with steady heat flux applied to one side and adiabatic and no-slip conditions applied to the opposite side, top and bottom.  The heated side is defined as a free surface, and is subject to radiative and convective losses. Gasification is not considered in this first example. Material properties except for viscosity are taken as constant, with values as given in Table I.
583
584
For polymer melts, viscosity is a function of both temperature and molecular weight distribution, which changes as the polymer bonds break during heating. This model uses an approach that incorporates this degradation into a simple description of viscosity as a function of temperature alone <span id='citeF-22'></span>[[#cite-22|[22]]]. Figure [[#img-5|5]] shows three curves of viscosity vs. temperature for the polypropylene type PP702N, a low viscosity commercial injection molding resin formulation.  The relationship used in the model, as shown by the black line, connects the curve for the undegraded polymer to points A and B extrapolated from the viscosity curve for each melt sample to the temperature at which the sample was formed.  The viscosity at room temperature is set to <math display="inline">10^{6}</math> Pa-s to maintain rigidity, and viscosity changes linearly between <math display="inline">T=25</math> <math display="inline">^\circ </math>C and <math display="inline">T=200</math> <math display="inline">^\circ </math>C. The result is an empirical viscosity-temperature curve that implicitly accounts for molecular weight changes. The analytical definition of this curve is given in Box II.
585
586
An initial spacing of 2.0 mm between nodes results in a finite element model of 1537 nodes and 2818 three-noded triangular elements for a 2D representation of the geometry, while a spacing of 1.4 mm results in 3098 nodes and 5832 elements. No nodes are added during the course of the analysis in this and all the other examples shown in the paper.
587
588
The addition of a catch pan to capture the dripping polymer melt tests the ability of the PFEM model to recover mass when a particle or set of particles reaches the catch surface.  For this problem, heat flux was only applied to free surfaces above the midpoint between the catch pan and the base of the sample. However, every free surface is subject to radiative and convective heat losses. The problem was studied for three heat flux levels of 20, 30 and 40 kW/m<math display="inline">^{2}</math>.  The catch pan is set to 250 <math display="inline">^\circ </math>C, a temperature that is high enough to keep the melt fluid.  The thickness of the sample is reduced to 2.5 cm to achieve results more quickly, and the height of 25 cm reflects the dimensions of an earlier experiment <span id='citeF-22'></span>[[#cite-22|[22]]].
589
590
Figure [[#img-7|7]]a shows four snapshots of the time evolution of the melt flow into the catch pan using the finer grid and a heat flux level of 20 kW/m<math display="inline">^{2}</math>.
591
592
<div id='img-7'></div>
593
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
594
|-
595
|[[Image:Draft_Samper_846581836-Fig4_polymers.png|600px|(a) Evolution of the melt flow into the catch pan at t = 400s, 550s, 700s and 1000s]]
596
|- style="text-align: center; font-size: 75%;"
597
| colspan="1" | '''Figure 7:''' (a) Evolution of the melt flow into the catch pan at t = 400s, 550s, 700s and 1000s
598
|}
599
600
<div id='img-8'></div>
601
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
602
|-
603
|[[Image:Draft_Samper_846581836-Fig4b_polymers1.png|312px|]]
604
|[[Image:Draft_Samper_846581836-Fig4b_polymers2.png|324px|(b) Detail of the flow of the melt as it drips into the pan at t=550 and 1000s]]
605
|- style="text-align: center; font-size: 75%;"
606
| colspan="2" | '''Figure 8:''' (b) Detail of the flow of the melt as it drips into the pan at <math>t=550</math> and 1000s
607
|}
608
609
A detail of the flow of the melt as it drips into the pan at two different times is shown in Figure [[#img-8|8]]b.
610
611
Figure [[#img-9|9]] (left) shows the flow into the catch pan at time <math display="inline">t = 600</math> s. Directly below the base of the sample where the melt is dripping, the temperature of the melt is higher. On the catch pan away from this point, the top of the melt has cooled below 325 <math display="inline">^\circ </math>C, the temperature of the catch pan surface in this case. The melt spreads to either side from the point at which the dripping melt contacts the catch pan.
612
613
Figure [[#img-9|9]] (right) shows the mass of the sample, the mass of the melt on the catch pan, and their sum.  After a heating time of about 170 s, the mass begins to be transferred from the sample to the catch pan.  The total mass reflects a conservation of mass within <math display="inline">\pm 5%</math>.  Note that because of the way nodes are eliminated and recreated at surfaces in the PFEM, the total mass at times exceeds 100% of the initial mass. Holes in the melt flow on the catch plate are occasionally generated by the Alpha-Shape method. Both types of error can be reduced by increasing the number of particles (nodes), as in the standard FEM.
614
615
<div id='img-9'></div>
616
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
617
|-
618
|[[Image:Draft_Samper_846581836-Fig5_polymers.png|600px|Melt flow into catch pan at t = 600s. Mass vs time for polymer in sample,   in catch pan, and total mass ]]
619
|- style="text-align: center; font-size: 75%;"
620
| colspan="1" | '''Figure 9:''' Melt flow into catch pan at t = 600s. Mass vs time for polymer in sample,   in catch pan, and total mass 
621
|}
622
623
Figure [[#img-10|10]] compares the mass loss rate of the quasi-steady period with that obtained from experiments at three levels of heat flux performed at NIST <span id='citeF-22'></span>[[#cite-22|[22]]]. The mass loss rate follows the same trend, although the values are about 25% higher than the experimental data. Note that this solution includes only melt flow in response to heating at a steady heat flux with radiative and convective losses.
624
625
===4.3 Inclusion of gasification effects===
626
627
The same problem was solved next ''including gasification effects''. The parameters in the gasification heat flux model chosen for the computations (Eqs.(13)&#8211;(15)) are given in Table I.
628
629
The general trend of the numerical results including gasification is not very different from the results shown in Figures [[#img-7|7]] and [[#img-9|9]]. Gasification effects however have an influence on the mass loss rate.
630
631
Results for the mass loss rate including gasification are also plotted in Figure [[#img-10|10]]. The computational results fall within experimental error for the range of heat fluxes considered. This shows the importance of accounting for gasification effects in this type of problem. The decrease in mass loss rate when gasification is added to flow demonstrates two competing effects:  the loss of additional mass as gas is released and the cooling caused by the endothermic process of polymer decomposition.
632
633
Note that this model is not yet complete, since phenomena such as in-depth absorption, the melting of the crystalline fraction of the polymer, and the temperature dependence of other material properties are not yet included. However, these phenomena are not difficult to add to the model and will be included in future development.
634
635
The solution of this melt flow problem took 2 CPU hours in a standard Pentium 4 PC.
636
637
<div id='img-10'></div>
638
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
639
|-
640
|[[Image:Draft_Samper_846581836-grafico.png|540px|Comparison of PFEM results to experiments for mass loss rate as a function of incident radiant flux]]
641
|- style="text-align: center; font-size: 75%;"
642
| colspan="1" | '''Figure 10:''' Comparison of PFEM results to experiments for mass loss rate as a function of incident radiant flux
643
|}
644
645
===4.4 3D polymer melt flow===
646
647
To test the ability of the PFEM to solve polymer melt flow  problems in three dimensions, a 3D heated sample was studied. The same boundary conditions are used as in the 2D problem illustrated in Figure [[#img-5|5]], but the initial dimensions of the sample are reduced to <math display="inline">10</math> cm <math display="inline">\times 10</math> cm <math display="inline">\times 2.5</math> cm thick.  The initial discretization has 22475 nodes and 97600 four-noded tetrahedra, and the runtime for this problem was slightly over 8 hours in a Pentium 4 PC.  The shape of the surface and temperature field at different times after heating begins are shown in Figure [[#img-11|11]].
648
649
<div id='img-11'></div>
650
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
651
|-
652
|[[Image:Draft_Samper_846581836-Fig7_polymers.png|527px|Solution of a 3D polymer melt problem with the PFEM. Melt flow from a   heated prismatic sample at different times.]]
653
|- style="text-align: center; font-size: 75%;"
654
| colspan="1" | '''Figure 11:''' Solution of a 3D polymer melt problem with the PFEM. Melt flow from a   heated prismatic sample at different times.
655
|}
656
657
Edge effects in the 3D model slow the rate of flow along the side walls, resulting in a thicker sample there throughout the melt flow process.  This has also been observed in the experiments.  Although the resolution for this problem is not fine enough to achieve high accuracy, the qualitative agreement of the 3D model with 2D flow and the ability to carry out this problem in a reasonable amount of time indicate that the PFEM can be used to model melt flow and spread of complex 3D polymer geometry.
658
659
===4.5 Melting and flow of a triangular slab===
660
661
Figure [[#img-12|12]] shows results for the analysis of the melting of a triangular thermoplastic object with a base of <math display="inline">5</math> cm and a height of <math display="inline">5\times 5</math> cms.  The solution is symmetrical about the <math display="inline">x=0</math> axis. All free surfaces are heated at a heat flux of 30 kW/m<math display="inline">^2</math>. The material properties for the polymer are the same as for the previous 2D example. The potential of the PFEM for simulating the progressive detachment and flow of the polymer particles from the object surface towards the underlying floor is clearly demonstrated.
662
663
A similar but more complex analysis is shown in Figure [[#img-13|13]]. Here the problem is the simulation of the melt flow of the same triangular thermoplastic object into a catch pan. The PFEM succeeds in predicting in a very realistic manner the progressive melting and slip of the polymer particles along the vertical wall separating the triangular object and the catch pan. The analysis follows until the whole object has fully melted and its mass is transferred to the catch pan. The total mass was preserved with an accuracy of 0.5% in both these studies. Gasification, in-depth absorption and radiation were not taken into account in these examples.
664
665
<div id='img-12'></div>
666
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
667
|-
668
|[[Image:Draft_Samper_846581836-Fig8_polymers.png|600px|Simulation of the melt flow of a heated triangular thermoplastic object]]
669
|- style="text-align: center; font-size: 75%;"
670
| colspan="1" | '''Figure 12:''' Simulation of the melt flow of a heated triangular thermoplastic object
671
|}
672
673
<div id='img-13'></div>
674
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
675
|-
676
|[[Image:Draft_Samper_846581836-Fig9_polymers.png|438px|Melt flow of a heated triangular object into a catch pan]]
677
|- style="text-align: center; font-size: 75%;"
678
| colspan="1" | '''Figure 13:''' Melt flow of a heated triangular object into a catch pan
679
|}
680
681
===4.6 Study of melt spread===
682
683
In the 2D PFEM model of the melt spread experiment, the initial computational space is an upright rectangle representing the thermoplastic object mounted above a rectangular catch plate, as illustrated in Figure [[#img-14|14]].  The catch plate may be tilted with respect to the horizontal line, although the catch plate always extends along the x-axis for ease of defining the distinction between the melt belonging to the object and that in the melt pool.
684
685
The left surface of the object and the top surface of the catch plate are designated as free surfaces, which are subject to heat losses from radiation and convection.  A steady heat flux is applied only to the free surface of the object, whose drips are distinguished from the melt on the catch plate by location above a specified height.  All other faces in the problem, identified by dark lines in Figure [[#img-14|14]], obey no-slip conditions. These faces are adiabatic, except for the bottom surface of the catch plate, which is maintained at a fixed heater temperature of <math display="inline">235</math> <math display="inline">^\circ </math>C. The thermoplastic object is initially at room temperature, and the linear temperature distribution within the catch plate balances the heat losses from the upper surface.  The temperature of the top surface, given the heater temperature and the thermal conductivity of the catch plate, was found to agree with the experimental value <span id='citeF-21'></span>[[#cite-21|[21]]].
686
687
Thermophysical material properties for the polymer are the same as for the other examples in this paper. The viscosity of the dripping object follows the curve for the undegraded polymer in Figure [[#img-5|5]]. For the flow of melt on the catch plate, a straight line is fitted to the viscosity-temperature relationship for this polymer after it has been degraded at a heat flux of 22 kW m<math display="inline">^2</math>. This relationship is obtained by rheological measurements of the melt deposited on the catch plate, and is not shown in Figure [[#img-5|5]]. The ceramic catch plate is assigned the thermal properties reported by the manufacturer. Neither gasification, in-depth absorption, nor radiation were taken into account in this study.
688
689
The spatial resolution of the finite element mesh is initially uniform throughout the thermoplastic object. Results are reported for a  grid with initial spacing of 1.0 mm.  For the catch plate, the mesh size at the top surface must be small enough to catch the particles dripping from the object, but the mesh can be coarser further in-depth since only thermal diffusion needs to be resolved.
690
691
<div id='img-14'></div>
692
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
693
|-
694
|[[Image:Draft_Samper_846581836-Fig10_polymers.png|529px|Initial conditions for melt spread problem]]
695
|- style="text-align: center; font-size: 75%;"
696
| colspan="1" | '''Figure 14:''' Initial conditions for melt spread problem
697
|}
698
699
Figure [[#img-15|15]] shows a snapshot of flow onto a horizontal catch plate at time <math display="inline">t = 600</math> s.  Directly below the base of the object where the polymer melt is dripping, the temperature of the melt pool is around 340 <math display="inline">^\circ </math>C. This equals the set temperature of the melt feeder used in the experiments. The temperature of the melt pool cools as it flows away from the line of dripping and loses heat to the catch plate and to the surroundings.  At far enough distances the surface temperature of the melt pool drops below the temperature of the catch plate.  If the plate was not heated, the polymer would continue to cool, until its viscosity may be high enough to interfere with the flow.
700
701
<div id='img-15'></div>
702
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
703
|-
704
|[[Image:Draft_Samper_846581836-Fig11_polymers.png|600px|Thermoplastic object dripping into horizontal catch pan at t = 600 s]]
705
|- style="text-align: center; font-size: 75%;"
706
| colspan="1" | '''Figure 15:''' Thermoplastic object dripping into horizontal catch pan at <math>t = 600</math> s
707
|}
708
709
Figure [[#img-16|16]] shows the evolution of the melt pool with time as the object continues to drip.  The spread rate is determined by a balance between the potential energy caused by the accumulation of melt under the drip line, the kinetic energy of the flow and the viscous drag on the melt as it moves over the catch plate.
710
711
<div id='img-16'></div>
712
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
713
|-
714
|[[Image:Draft_Samper_846581836-Fig12_polymers.png|600px|Spread of melt pool over horizontal catch plate]]
715
|- style="text-align: center; font-size: 75%;"
716
| colspan="1" | '''Figure 16:''' Spread of melt pool over horizontal catch plate
717
|}
718
719
In Figure [[#img-17|17]] the catch plate has been tilted at an angle of 1.8 degrees.  The slope of the plate adds to the potential energy of the melt, and the spread rate increases in the downslope direction while it decreases in the upstream direction where the accumulation of fluid must counter the potential energy due to the slope.  A study of the movement of the melt front with time for horizontal and tilted catch plates is reported in <span id='citeF-21'></span>[[#cite-21|[21]]]. PFEM spread rates are within <math display="inline">10 %</math> of experimental results when the continuing degradation over the surface of the heated catch plate is taken into account.
720
721
<div id='img-17'></div>
722
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
723
|-
724
|[[Image:Draft_Samper_846581836-Fig13_polymers.png|600px|Thermoplastic object dripping into tilted catch plate at t = 600 s]]
725
|- style="text-align: center; font-size: 75%;"
726
| colspan="1" | '''Figure 17:''' Thermoplastic object dripping into tilted catch plate at <math>t = 600</math> s
727
|}
728
729
===4.7 Melting of a thermoplastic chair===
730
731
This example and the example that follows demonstrate the capabilities of the PFEM to model the fire behavior of polymers of more complex shapes.
732
733
The example shown in Figure [[#img-18|18]] is the simulation of the melt flow of a thermoplastic object resembling a chair modeled as a 2D solid. The images show the progressive melting of the chair exposed to a heat flux of 30 kW/m<math display="inline">^2</math> on all surfaces except the base. The ability of the PFEM to model self-contact situations as the shape of the chair changes with time due to melting is demonstrated.
734
735
<div id='img-18'></div>
736
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
737
|-
738
|[[Image:Draft_Samper_846581836-Fig14_polymers.png|600px|Melting of a heated chair modeled as a 2D object. Note self-contact between chair surfaces as melting evolves]]
739
|- style="text-align: center; font-size: 75%;"
740
| colspan="1" | '''Figure 18:''' Melting of a heated chair modeled as a 2D object. Note self-contact between chair surfaces as melting evolves
741
|}
742
743
<div id='img-19'></div>
744
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
745
|-
746
|[[Image:Draft_Samper_846581836-Fig15_polymers.png|600px|Melting of a cylindrical candle with the PFEM]]
747
|- style="text-align: center; font-size: 75%;"
748
| colspan="1" | '''Figure 19:''' Melting of a cylindrical candle with the PFEM
749
|}
750
751
===4.8 Melting of a candle===
752
753
The last example shown in Figure [[#img-19|19]] and [[#img-20|20]] is the simulation of the melting of a cylindrical candle. Heat flux is applied at the top surface of the candle at a rate of 30 KW/m<math display="inline">^{2}</math>. This induces the progressive melting and dripping of the candle particles along the cylindrical wall until they eventually hit and spread over the underlying pan surface.
754
755
<div id='img-20'></div>
756
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
757
|-
758
|[[Image:Draft_Samper_846581836-Fig16_polymers.png|450px|Melting of candle. Detail of melted zone and pan surface at a certain instant ]]
759
|- style="text-align: center; font-size: 75%;"
760
| colspan="1" | '''Figure 20:''' Melting of candle. Detail of melted zone and pan surface at a certain instant 
761
|}
762
763
==5 CONCLUDING REMARKS==
764
765
The PFEM is a powerful technique to model the melting, flow and flame spread of thermoplastic objects in fire situations. The method allows to track the motion of the polymer particles as they melt, flow over the surface of the object and fall towards and on the underlying floor. The PFEM can also predict the spread of the flame in the floor for different ambient temperature conditions and effect of gasification, in-depth absorption and radiation problems.
766
767
The simulation of the melting of a simplified chair of arbitrary shape and a candle has shown the potential of the PFEM to model the drastic change of shape of objects as they melt, drip and spread, including self-contact situations.
768
769
Developments in progress include coupling the PFEM formulation with the surrounding hot air and the external heat source, inclusion of elastic effects in the polymer material and increasing the computational efficiency of the method.
770
771
==ACKNOWLEDGMENTS==
772
773
The authors thank Dr. T. Ohlemiller at NIST for providing experimental results for the 2D rectangular slab problem <span id='citeF-22'></span>[[#cite-22|[22]]]. Special thanks are given to Mr. P. Ryzhakov and to Drs. J. Marti and F. del Pin from CIMNE for their advice during this research. The support of Dr. P. Dadvand in the development of the PFIRE code using the Kratos software platform <span id='citeF-23'></span>[[#cite-23|[23]]] is acknowledged.
774
775
===BIBLIOGRAPHY===
776
777
<div id="cite-1"></div>
778
'''[[#citeF-1|[1]]]'''  Zhang J, Shields TJ, Silcock GWH.  Effect of melting behaviour on flame spread of thermoplastics. ''Fire & Materials'' 1997; '''21'''(1):1&#8211;6.
779
780
<div id="cite-2"></div>
781
'''[2]'''  Fleischmann CM, Hill GR. Burning behaviour of upholstered furniture. Interflam 2004; 907&#8211;916.
782
783
<div id="cite-3"></div>
784
'''[[#citeF-3|[3]]]'''  Sherratt J, Drysdale D. The effect of the melt-flow process on the fire behaviour of thermoplastics. ''Interflam'' 2001; 149&#8211;159.
785
786
<div id="cite-4"></div>
787
'''[[#citeF-4|[4]]]'''  Butler KM, Ohlemiller TJ, Linteris GT. A progress report on numerical modeling of experimental polymer melt flow behavior. ''Interflam'' 2004; pp. 937&#8211;948.
788
789
<div id="cite-5"></div>
790
'''[[#citeF-5|[5]]]'''  Idelsohn SR, Oñate E, Del Pin F.  A lagrangian meshless finite element method applied to fluid-structure interaction problems. ''Computer and Structures'' 2003b; '''81''':655&#8211;671.
791
792
<div id="cite-6"></div>
793
'''[[#citeF-6|[6]]]'''  Idelsohn SR, Oñate E, Del Pin F.   The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. ''International Journal for Numerical  Methods in Engineering'' 2004; '''61''':964&#8211;989.
794
795
<div id="cite-7"></div>
796
'''[7]'''  Idelsohn SR, Oñate E, Calvo N, Del Pin F.   The meshless finite element method.  ''International Journal for Numerical  Methods in Engineering'' 2003a; '''58'''(6): 893&#8211;912.
797
798
<div id="cite-8"></div>
799
'''[[#citeF-8|[8]]]'''  Oñate E, Idelsohn SR, Del Pin F, Aubry R.  The particle finite element method. An overview. ''International Journal of Computational Methods'' 2004; '''1'''(2):267&#8211;307.
800
801
<div id="cite-9"></div>
802
'''[[#citeF-9|[9]]]'''  Aubry R, Idelsohn SR, Oñate E.  Particle finite element method in fluid-mechanics including thermal convection-diffusion. ''Computers and Structures'' 2005; '''83'''(17-18):1459&#8211;1475.
803
804
<div id="cite-10"></div>
805
'''[[#citeF-10|[10]]]'''  Aubry R, Idelsohn SR, Oñate E.   Fractional step like schemes for free surface problems with thermal coupling using the Lagrangian PFEM. ''Computational Mechanics'' 2006; '''38'''(4-5):294&#8211;309.
806
807
<div id="cite-11"></div>
808
'''[[#citeF-11|[11]]]'''    Oñate E, Idelsohn SR, Celigueta MA, Rossi R.  Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. ''Computer Methods in Applied Mechanics and Engineering'' 2008; '''197'''(19-20):1777&#8211;1800.
809
810
<div id="cite-12"></div>
811
'''[[#citeF-12|[12]]]'''  Butler KM, Oñate E, Idelsohn SR, Rossi R. Modeling polymer melt flow using the particle finite element method. 11th International Interflam Conference (''Interflam'07''), September 3&#8211;5 2007, London, England, Volume 2,  pp. 929&#8211;940.
812
813
<div id="cite-13"></div>
814
'''[[#citeF-13|[13]]]'''  Rossi R, Butler KM, Oñate E,  Idelsohn SR. Modeling polymer melt flow using the Particle Finite Element Method. ''Advanced Research Workshop on Fire Computer Modeling'', Gijón, Spain, October 18-20, 2007.
815
816
<div id="cite-14"></div>
817
'''[[#citeF-14|[14]]]'''  Edelsbrunner H, Mucke EP.  Three dimensional alpha shapes. ''ACM Trans. Graphics'' 1999; '''13''':43&#8211;72.
818
819
<div id="cite-15"></div>
820
'''[[#citeF-15|[15]]]'''  Zienkiewicz OC, Taylor RL, Nithiarasu N. ''The Finite Element Method. Vol. 3 Fluid Mechanics'', Elsevier, 2005.
821
822
<div id="cite-16"></div>
823
'''[[#citeF-16|[16]]]'''  Oñate E. Derivation of stabilized equations for numerical solution of advective-diffusive transport and fluid flow problems. ''Computer Methods in Applied Mechanics and Engineering'' 1998; '''151''':233-265.
824
825
<div id="cite-17"></div>
826
'''[[#citeF-17|[17]]]'''  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Computer Methods in Applied Mechanics and Engineering'' 2000; '''182'''(3-4):355-370.
827
828
<div id="cite-18"></div>
829
'''[[#citeF-18|[18]]]'''  Oñate E. Possibilities of finite calculus in computational mechanics. ''International Journal for Numerical Methods in Engineering'' 2004; '''60'''(1):255&#8211;281.
830
831
<div id="cite-19"></div>
832
'''[[#citeF-19|[19]]]'''  Zienkiewicz OC, Jain P C, Oñate E.  Flow of solids during forming and extrusion: Some aspects of numerical solutions. ''International Journal of Solids and Structures'' 1978; '''14''':15&#8211;38.
833
834
<div id="cite-20"></div>
835
'''[[#citeF-20|[20]]]'''  Zienkiewicz OC, Oñate E, Heinrich JC.  A general formulation for the coupled thermal flow of metals using finite elements. ''International Journal for Numerical Methods in Engineering'' 1981; '''17''':1497&#8211;1514.
836
837
<div id="cite-21"></div>
838
'''[[#citeF-21|[21]]]'''  Butler KM.   A model of melting and dripping thermoplastic objects in fire. Fire and Materials 2009, San Francisco, USA, January 26&#8211;28, 2009, pp. 341&#8211;352.
839
840
<div id="cite-22"></div>
841
'''[[#citeF-22|[22]]]'''  Ohlemiller TJ, Shields JR.  Aspects of the Fire Behavior of Thermoplastic Materials.  ''NIST Technical Note 1493'' 2008.
842
843
<div id="cite-23"></div>
844
'''[[#citeF-23|[23]]]'''  Davdand P, Rossi R, Oñate E. An object-oriented environment for developing finite element codes for multi-disciplinary applications. ''Archives of Computational Methods in Engineering'' 2009, accepted for publication.
845
846
==APPENDIX==
847
848
The  matrices and vectors in Eqs.(9)-(12) for a 4-noded tetrahedron are:
849
850
{| class="formulaSCP" style="width: 100%; text-align: left;" 
851
|-
852
| 
853
{| style="text-align: left; margin:auto;width: 100%;" 
854
|-
855
| style="text-align: center;" | <math>\mathbf{M}_{ij} =\int _{V^{e}} \rho \mathbf{N}_{i}^{T} \mathbf{N}_{j} dV \quad , \quad \mathbf{K}_{ij}=\int _{V^{e}} \mathbf{B}_{i}^{T} \mathbf{D} \mathbf{B}_{j} dV </math>
856
|}
857
|}
858
859
{| class="formulaSCP" style="width: 100%; text-align: left;" 
860
|-
861
| 
862
{| style="text-align: left; margin:auto;width: 100%;" 
863
|-
864
| style="text-align: center;" | <math> \mathbf{G}_{ij} = \int _{V^{e}} \mathbf{B}_{i}^{T}  \mathbf{m}\mathbf{N}_{j} dV \quad , \quad \mathbf{f}_{i}=\int _{V^{e}} \mathbf{ N}_{i}^{T}\mathbf{b}dV+ \int _{\Gamma ^e} \mathbf{ N}_{i}^{T}\mathbf{t}d\Gamma </math>
865
|}
866
|}
867
868
{| class="formulaSCP" style="width: 100%; text-align: left;" 
869
|-
870
| 
871
{| style="text-align: left; margin:auto;width: 100%;" 
872
|-
873
| style="text-align: center;" | <math>{L}_{ij} = \int _{V^e} {\boldsymbol \nabla }^T N_i \tau {\boldsymbol \nabla }N_j dV  \quad ,\quad {\boldsymbol \nabla } = \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2},{\partial  \over \partial x_3}\right]^T</math>
874
|}
875
|}
876
877
{| class="formulaSCP" style="width: 100%; text-align: left;" 
878
|-
879
| 
880
{| style="text-align: left; margin:auto;width: 100%;" 
881
|-
882
| style="text-align: center;" | <math>{\boldsymbol Q}= [{\boldsymbol Q}_1,{\boldsymbol Q}_2,{\boldsymbol Q}_3]\quad ,\quad [\mathbf{Q}_{k}]_{ij} = \int _{V^e}\tau _k {\partial N_i \over \partial x_k} \mathbf{N}_j dV\quad ,\quad \hbox{no sum in }k</math>
883
|}
884
|}
885
886
{| class="formulaSCP" style="width: 100%; text-align: left;" 
887
|-
888
| 
889
{| style="text-align: left; margin:auto;width: 100%;" 
890
|-
891
| style="text-align: center;" | <math>\hat {\boldsymbol M}= \left[\hat {\boldsymbol M}_1, \hat {\boldsymbol M}_2,\hat {\boldsymbol M}_3\right]\,,\, [\hat {\boldsymbol M}_i]_{kl}=\left(\int _{V^e} \tau _k {\boldsymbol N}_i{\boldsymbol N}_j dV \right)\delta _{kl} \quad , \quad k,l=1,2,3</math>
892
|}
893
|}
894
895
{| class="formulaSCP" style="width: 100%; text-align: left;" 
896
|-
897
| 
898
{| style="text-align: left; margin:auto;width: 100%;" 
899
|-
900
| style="text-align: center;" | <math>{\boldsymbol B}=[{\boldsymbol B}_1,{\boldsymbol B}_2,{\boldsymbol B}_3,{\boldsymbol B}_4];\, {\boldsymbol B}_i =\left[\begin{matrix} \displaystyle {\partial N_i \over \partial x}&0&0\\ 0 & \displaystyle {\partial N_i \over \partial y} &0\\ 0&0& \displaystyle {\partial N_i \over \partial z} \\ \displaystyle {\partial N_i \over \partial y} & \displaystyle {\partial N_i \over \partial x}&0\\\\ \displaystyle {\partial N_i \over \partial z} &0&\displaystyle {\partial N_i \over \partial x}\\\\ 0& \displaystyle {\partial N_i \over \partial z} &\displaystyle {\partial N_i \over \partial y} \end{matrix} \right];\,\, {\boldsymbol D} =\mu \left[\begin{matrix}2 &0&0&0&0&0\\ 0&2&0&0&0&0\\ 0&0&2&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1\end{matrix}\right]</math>
901
|}
902
|}
903
904
{| class="formulaSCP" style="width: 100%; text-align: left;" 
905
|-
906
| 
907
{| style="text-align: left; margin:auto;width: 100%;" 
908
|-
909
| style="text-align: center;" | <math>{\boldsymbol N}=[{\boldsymbol N}_1,{\boldsymbol N}_2,{\boldsymbol N}_3,{\boldsymbol N}_4] \quad , \quad  {\boldsymbol N}_i={N}_i {\boldsymbol I}_3\quad ,\quad  {\boldsymbol I}_3 \quad \hbox{is the }3\times 3 \hbox{ unit matrix}</math>
910
|}
911
|}
912
913
{| class="formulaSCP" style="width: 100%; text-align: left;" 
914
|-
915
| 
916
{| style="text-align: left; margin:auto;width: 100%;" 
917
|-
918
| style="text-align: center;" | <math>{C}_{ij}=\int _{V^e} \rho {c}{N}_i {N}_j dV \quad ,\quad  {H}_{ij}= \int _{V^e}{\boldsymbol \nabla }^T {\boldsymbol N}_i [k] {\boldsymbol \nabla }{\boldsymbol N}_j dV</math>
919
|}
920
|}
921
922
{| class="formulaSCP" style="width: 100%; text-align: left;" 
923
|-
924
| 
925
{| style="text-align: left; margin:auto;width: 100%;" 
926
|-
927
| style="text-align: center;" | <math>[k]=  \left[\begin{matrix}k_1 &0 &0 \\ 0 & k_2&0\\ 0&0& k_3\end{matrix}\right]\quad , \quad   q_i =\int _{V^e}N_i (Q+q_{gas}) dV - \int _{\Gamma _q^{e}}  N_i \bar q_n d\Gamma </math>
928
|}
929
|}
930
931
In above equations indexes <math display="inline">i,j</math> run from 1 to the number of element nodes (4 for a tetrahedron), <math display="inline">q_n</math> is the heat flow prescribed at the external boundary <math display="inline">\Gamma </math>, '''t''' is the surface traction vector <math display="inline">\mathbf{t}=[t_x,t_y,t_z]^T</math> and <math display="inline">V^{e}</math> and <math display="inline">\Gamma ^e</math> are the element volume and the element boundary, respectively.
932

Return to Onate et al 2010a.

Back to Top