Melting and Spread of Polymers in Fire with the Particle Finite Element Method

Published in ''Int. Journal for Numerical Methods in Engineering'' Vol. 81 (8), pp. 1046-1072, 2010<br />
doi: 10.1002/nme.2731

==Abstract==

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.

'''Keywords''': Melting, dripping, polymers, Particle Finite Element Method, PFEM

==1 INTRODUCTION==

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.

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]]].

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.

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]]].

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.

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]]].

==2 THE PARTICLE FINITE ELEMENT METHOD. AN OVERVIEW==

===2.1 Basic steps of the PFEM===

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.

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.

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.

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>).

A typical solution with the PFEM involves the following steps:

<ol>

<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>
<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>
<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>
<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>
<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>
<li>Go back to step 1 and repeat the solution process for the next time step. </li>

</ol>

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Figure1.png|500px|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]]
|- style="text-align: center; font-size: 75%;"
| 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
|}

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.

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]]].

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]]].

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig2a.png|540px|]]
|[[Image:Draft_Samper_846581836-Fig2b.png|540px|Modelling of contact between the melting object and a fixed boundary]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Modelling of contact between the melting object and a fixed boundary
|}

===2.2 Governing equations===

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:

==Momentum==

<span id="eq-1"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho \frac{dv_i}{dt}= \frac{\partial \sigma _{ij}}{\partial x_j}+b_i \quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

==Mass balance==

<span id="eq-2"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial v_{i}}{\partial x_i}=0 \quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

==Heat transport==

<span id="eq-3"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

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)

<span id="eq-4"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij} =s_{ij}-p\delta _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4.a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4.b)
|}

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>.

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.

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).

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:

<span id="eq-5.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k\frac{\partial T}{\partial n}+\bar q_n =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.a)
|}

with

<span id="eq-5.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.b)
|}

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.

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.

===2.3 Discretization of the equations===

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

<span id="eq-6"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">\tau </math> is a stabilization parameter given by <span id='citeF-17'></span>[[#cite-17|[17]]]

<span id="eq-7"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = \left(\frac{2\rho \vert \mathbf{v}\vert }{h}+\frac{8\mu }{3h^2} \right)^{-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

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

<span id="eq-8"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">w_i</math> are arbitrary weighting functions.

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]]].

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

==Momentum==

<span id="eq-9"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

==Mass balance==

<span id="eq-10"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

==Pressure gradient projections==

<span id="eq-11"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat {\boldsymbol M} \bar {\boldsymbol \pi } + {\boldsymbol Q}^T \bar {\boldsymbol p} ={\boldsymbol 0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

==Heat transfer==

<span id="eq-12"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol C}\dot{\bar{\boldsymbol T}} +{\boldsymbol H}\bar {\boldsymbol T}= {\boldsymbol q} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

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.

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.
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]]].

==3 ACCOUNTING FOR GASIFICATION EFFECTS==

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

1. LOOP OVER TIME STEPS, <math>t=1</math>, NTIME

Known values

<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>


2. LOOP OVER NUMBER OF ITERATIONS, <math>i=1</math>, NITER

<math>\bullet </math> Compute the nodal velocities by solving Eq. [[#eq-9|(9)]]

<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>

<math>\bullet </math> Compute  nodal pressures from Eq.[[#eq-10|(10)]]

<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>

<math>\bullet </math> Compute nodal pressure gradient projections from Eq. [[#eq-11|(11)]]

<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>

<math>\bullet </math> Compute nodal temperatures from Eq.[[#eq-12|(12)]]

<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>

<math>\bullet </math> Update position of analysis domain nodes:

<math>{~}^{t+1}\bar {\boldsymbol x}^{i+1}= {~}^t{\boldsymbol x}^i+ {~}^{t+1}{\boldsymbol  v}^{i+1}\Delta t</math>


Define new “cloud” of nodes <math>^{t+1} {C}^{i+1}</math>

<math>\bullet </math> Update viscosity values in terms of temperature

<math>{~}^{t+1}\mu =\mu (^{t+1}\bar {\boldsymbol T}^{i+1})</math>


Check convergence  for all variables <math>\to </math> NO <math>\to </math>  Next iteration <math>i\to i+1</math>

<math>~~~~~~~~~~\downarrow </math> YES


Next time step <math>t\to t+1</math>

<math>\bullet </math> Identify new analysis domain boundary: <math>^{t+1}V</math>

<math>\bullet </math> Generate mesh: <math>{~}^{t+1}M</math>

Go to 1

'''Box I'''. Flow chart of basic PFEM algorithm for the polymer domain

<span id="eq-13"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{gas}=\rho H\bar \varepsilon _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">H</math> is the enthalpy of vaporization,

<span id="eq-14"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar \varepsilon _v =f(T) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

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]]]

<span id="eq-15"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(T) = - Ae^{-\frac{E}{RT}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

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.

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.

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.

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

<span id="eq-16"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

The new discretized mass balance equation is

<span id="eq-17"></span> 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

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.

An example of the effects of gasification will be shown in the following section.

==4 NUMERICAL EXAMPLES==

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]]].

===4.1 Validation of the gasification model===

A simple computational test  has been chosen to verify the accuracy of the gasification algorithm.

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 [[#table-1|I]]. The analysis was carried out with several meshes of 3-noded triangles.

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]]].

The results of  mass versus time for the different meshes considered are represented in Figure [[#img-3|3]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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.]]
|- style="text-align: center; font-size: 75%;"
| 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.
|}

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.

Some screenshots of the sample geometry at different times  for the finer mesh are shown in Figure [[#img-4|4]].

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.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-malla_gasification.png|360px|]]
|[[Image:Draft_Samper_846581836-temp_gasification_2000.png|420px|]]
|-
|[[Image:Draft_Samper_846581836-temp_gasification_4000.png|420px|]]
|[[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.]]
|- style="text-align: center; font-size: 75%;"
| 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.
|}

===4.2 Melting and flow of a rectangular slab===

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.

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 [[#cite-4|[4,?,13]]].
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.

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.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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.]]
|- style="text-align: center; font-size: 75%;"
| 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.
|}

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.

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]]].

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ 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
|-
| style="text-align: left;" |    Parameter 
| Symbol 
| Value 
| Units
|-
| style="text-align: left;" |    Density 
| <math>\rho </math>
| 900 
| kg/m<math display="inline">^3</math> 
|-
| style="text-align: left;" | Thermal conductivity 
| <math>k</math>
| 0.25 
| W/m-K  
|-
| style="text-align: left;" | Specific heat 
| <math>c</math>
| 2400 
| J/kg-K  
|-
| style="text-align: left;" | Reference temperature 
| <math>T_0</math>
| 298 
| K  
|-
| style="text-align: left;" | Emissivity 
| <math>\epsilon </math>
| 1.0 
| - 
|-
| style="text-align: left;" | Convective heat transfer coefficient 
| <math>\alpha _c</math>
| 8 
| W/m<math display="inline">^2</math>-K 
|-
| style="text-align: left;" | Heat of vaporization 
| <math>H</math>
| <math>8\times 10^5</math>
| J/kg  
|-
| style="text-align: left;" | Pre-exponential function 
| <math>A</math>
| <math>2.18\times 10^{12}</math>
| s<math display="inline">^{-1}</math> 
|-
| style="text-align: left;" | Activation energy divided by universal
|-
| style="text-align: left;" | gas constant 
| <math>E/R</math>
| 24400 
| K 
|-
| style="text-align: left;" | Stefan-Boltzmann constant 
| <math>\sigma </math>
| 5.67<math display="inline">\times 10^{-8}</math>
| W/m<math display="inline">^{2}</math>- K<math display="inline">^{4}</math>

|}

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-

<pre>

Tc: Temperature in degrees Celsius; T: Temperature in Kelvins

def f1(Tc):
return 10**(14.48 - 0.13858*Tc + 5.5960e-4*Tc*Tc - 7.8665e-7*Tc*Tc*Tc)

def f2(Tc):
return 10**(53.19 - 0.2542*Tc + 2.9879e-4*Tc*Tc)

def AuxFunction(T):

Tc = T - 273
if( Tc <= 25):
mu = 1e6
elif (Tc > 25 and  Tc <= 200):
mu = 1e6*(200-Tc)/175 + f1(200)
elif(Tc > 200 and Tc <= 350):
mu = f1(Tc)
elif(Tc > 350 and Tc < 425 ):
mu = f2(Tc)
else:
mu = f2(425)

return mu
</pre>


|}

Figure [[#img-6.a|6a]] 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>.
A detail of the flow of the melt as it drips into the pan at two different times is shown in Figure [[#img-6.b|6b]].

<div id='img-6.a'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6.(a):''' Evolution of the melt flow into the catch pan at t = 400s, 550s, 700s and 1000s
|}

<div id='img-6.b'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig4b_polymers1.png|312px|]]
|[[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]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6. (b):''' Detail of the flow of the melt as it drips into the pan at <math>t=550</math> and 1000s
|}

Figure [[#img-7|7]] (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.

Figure [[#img-7|7]] (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.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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 ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Melt flow into catch pan at t = 600s. Mass vs time for polymer in sample,   in catch pan, and total mass 
|}

Figure [[#img-8|8]] 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.

===4.3 Inclusion of gasification effects===

The same problem was solved next ''including gasification effects''. The parameters in the gasification heat flux model chosen for the computations (Eqs.([[#eq-13|13]]) &#8211; ([[#eq-15|15]]) are given in Table [[#table-1|I]].

The general trend of the numerical results including gasification is not very different from the results shown in Figures [[#img-6|6]] and [[#img-7|7]]. Gasification effects however have an influence on the mass loss rate.

Results for the mass loss rate including gasification are also plotted in Figure [[#img-8|8]]. 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.

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.

The solution of this melt flow problem took 2 CPU hours in a standard Pentium 4 PC.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-grafico.png|540px|Comparison of PFEM results to experiments for mass loss rate as a function of incident radiant flux]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Comparison of PFEM results to experiments for mass loss rate as a function of incident radiant flux
|}

===4.4 3D polymer melt flow===

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-9|9]].

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Solution of a 3D polymer melt problem with the PFEM. Melt flow from a   heated prismatic sample at different times.
|}

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.

===4.5 Melting and flow of a triangular slab===

Figure [[#img-10|10]] 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.

A similar but more complex analysis is shown in Figure [[#img-11|11]]. 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.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig8_polymers.png|600px|Simulation of the melt flow of a heated triangular thermoplastic object]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Simulation of the melt flow of a heated triangular thermoplastic object
|}

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig9_polymers.png|438px|Melt flow of a heated triangular object into a catch pan]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Melt flow of a heated triangular object into a catch pan
|}

===4.6 Study of melt spread===

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-12|12]].  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.

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-12|12]], 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]]].

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.

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.

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig10_polymers.png|529px|Initial conditions for melt spread problem]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Initial conditions for melt spread problem
|}

Figure [[#img-13|13]] 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.

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig11_polymers.png|600px|Thermoplastic object dripping into horizontal catch pan at t = 600 s]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Thermoplastic object dripping into horizontal catch pan at <math>t = 600</math> s
|}

Figure [[#img-14|14]] 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.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig12_polymers.png|600px|Spread of melt pool over horizontal catch plate]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Spread of melt pool over horizontal catch plate
|}

In Figure [[#img-15|15]] 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.

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig13_polymers.png|600px|Thermoplastic object dripping into tilted catch plate at t = 600 s]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Thermoplastic object dripping into tilted catch plate at <math>t = 600</math> s
|}

===4.7 Melting of a thermoplastic chair===

This example and the example that follows demonstrate the capabilities of the PFEM to model the fire behavior of polymers of more complex shapes.

The example shown in Figure [[#img-16|16]] 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.

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Melting of a heated chair modeled as a 2D object. Note self-contact between chair surfaces as melting evolves
|}

===4.8 Melting of a candle===

The last example shown in Figure [[#img-17|17]] and [[#img-18|18]] 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.

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig15_polymers.png|600px|Melting of a cylindrical candle with the PFEM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Melting of a cylindrical candle with the PFEM
|}

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_846581836-Fig16_polymers.png|450px|Melting of candle. Detail of melted zone and pan surface at a certain instant ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Melting of candle. Detail of melted zone and pan surface at a certain instant 
|}

==5 CONCLUDING REMARKS==

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.

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.

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.

==ACKNOWLEDGMENTS==

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.

===REFERENCES===

<div id="cite-1"></div>
'''[[#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.

<div id="cite-2"></div>
'''[2]'''  Fleischmann CM, Hill GR. Burning behaviour of upholstered furniture. Interflam 2004; 907&#8211;916.

<div id="cite-3"></div>
'''[[#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.

<div id="cite-4"></div>
'''[[#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.

<div id="cite-5"></div>
'''[[#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.

<div id="cite-6"></div>
'''[[#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.

<div id="cite-7"></div>
'''[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.

<div id="cite-8"></div>
'''[[#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.

<div id="cite-9"></div>
'''[[#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.

<div id="cite-10"></div>
'''[[#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.

<div id="cite-11"></div>
'''[[#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.

<div id="cite-12"></div>
'''[[#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.

<div id="cite-13"></div>
'''[[#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.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Edelsbrunner H, Mucke EP.  Three dimensional alpha shapes. ''ACM Trans. Graphics'' 1999; '''13''':43&#8211;72.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Zienkiewicz OC, Taylor RL, Nithiarasu N. ''The Finite Element Method. Vol. 3 Fluid Mechanics'', Elsevier, 2005.

<div id="cite-16"></div>
'''[[#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.

<div id="cite-17"></div>
'''[[#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.

<div id="cite-18"></div>
'''[[#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.

<div id="cite-19"></div>
'''[[#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.

<div id="cite-20"></div>
'''[[#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.

<div id="cite-21"></div>
'''[[#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.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  Ohlemiller TJ, Shields JR.  Aspects of the Fire Behavior of Thermoplastic Materials.  ''NIST Technical Note 1493'' 2008.

<div id="cite-23"></div>
'''[[#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.

==APPENDIX==

The  matrices and vectors in Eqs.([[#eq-9|9]]) &#8211; ([[#eq-12|12]]) for a 4-noded tetrahedron are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
|}

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.