A Computational Model for

Fire Growth and Spread

On Thermoplastic Objects

Michael Bockelie 

Qing Tang 

Reaction Engineering International

77 West 200 South

Suite 210

Salt Lake City, UT 84101 

NIST GCR 07-914

NIST GCR 07-914

A Computational Model for

Fire Growth and Spread

On Thermoplastic Objects

Prepared for 

U.S. Department of Commerce 

Building and Fire Research Laboratory 

National Institute of Standards and Technology 

Gaithersburg, MD 20899-8665

By 

Michael Bockelie 

Qing Tang 

Reaction Engineering International 

77 West 200 South

Suite 210

Salt Lake City, UT 84101

November 2007

U.S. Department of Commerce 

Carlos M. Gutierrez, Secretary 

National Institute of Standards and Technology 

 James M. Turner, Acting Director 
Notice

This report was prepared for the Building and Fire Research Laboratory

of the National Institute of Standards and Technology under Contract number

SB1341-05-C-0041. The statement and conclusions contained in this report

are those of the authors and do not necessarily reflect the views of the 

National Institute of Standards and Technology or the Building and Fire

Research Laboratory.
Committed Individuals Solving Challenging Problems

 

A Computational Model for Fire 

Growth and Spread on Thermoplastic Objects

Final Report 

By 

Mike Bockelie and Qing Tang 

Reaction Engineering International 

Salt Lake City, Utah 

October 31, 2007 

Prepared for 

National Institute of Standards and Technology 

Contract No. SB1341-05-C-0041 

REI Project No. 4716 

INTERNATIONAL

REACTION

ENGINEERING

77 West 200 South, Suite 210 - Salt Lake City, Utah 84101 

Telephone: (801) 364-6925 FAX: (801) 364-6977 

http://www.reaction-eng.com 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects i 

DISCLAMER 

“This report was prepared as an account of work sponsored by an agency of the United States 

Government. Neither the United States Government nor any agency thereof, nor any of their 

employees, makes any warranty, express or implied, or assumes any legal liability or 

responsibility for the accuracy, completeness, or usefulness of any information, apparatus, 

product, or process disclosed, or represents that its use would not infringe privately owned rights. 

Reference herein to any specific commercial product, process, or service by trade name, 

trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, 

recommendation, or favoring by the United States Government or any agency thereof. The views 

and opinions of authors expressed herein do not necessarily state or reflect those of the United 

States Government or any agency thereof.” 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects ii 

EXECUTIVE SUMMARY 

In this report is described the work effort by Reaction Engineering International (REI) to 

develop, demonstrate and deliver to the National Institute of Standards and Technology (NIST) a 

condensed phase computational fluid dynamics (CFD) based tool to model the processes of 

melting, flow and gasification of thermoplastic materials exposed to a high heat flux. Potential 

applications of the tool include investigating the behavior of polymer materials commonly used 

in personal computers and computer monitors if exposed to an intense heat flux, such as occurs 

during a fire 

The model delivered to NIST is based on a time dependent (time varying) grid CFD method. 

• The model is written in FORTRAN 90 in an object-oriented form. A 3D, finite volume, 

multi-block body-fitted time dependent (time varying) grid formulation is used to solve 

the unsteady Navier Stokes equations. The time integration, spatial discretization and 

overall solution procedure are based on standard CFD methods from the literature. A 

multi-grid method is used to accelerate convergence at each time step. 

• Sub-models are included to describe the temperature dependent viscosity relationship and 

in-depth gasification and absorption of thermoplastic materials, free surface flows and 

surface tension. NIST data is used for key material properties of the thermoplastic 

materials of interest. 

• A variety of boundary conditions can be used for the velocity field (no-slip, free-slip) and 

heat transfer to the object (adiabatic, heat loss, specified heat flux). 

• Model outputs include the time dependent velocity, temperature and position 

(displacement) at points in the thermoplastic body which can be imported to standard 

CFD visualization packages. Additional outputs include the time history of the mass loss 

rate and heat fluxes. 

• The accuracy and capabilities of the modeling tool are demonstrated on a series of test 

cases of increasing complexity. The test cases include grid sensitivity studies, adding heat 

loss boundary conditions, simulations for two thermoplastic materials (PP702N, PP6523), 

different heat flux scenarios and test problem configurations. 

Comparisons of model results to NIST experimental data indicate discrepancies between the 

model and experimental results, particularly for the rate of mass loss. To match NIST data for 

mass loss rate large changes were required to the parameters originally used in the in-depth 

absorption model and kinetic rate parameters in the in-depth gasification model. In addition, for 

simulations in which the free surface of the melt flow is parallel to the direction of gravity a grid 

smoothing operation needed to be applied to the free surface to control grid skewness that would 

lead to simulation divergence. For simulations that did not include models for all key physics 

(e.g., no in-depth absorption, no in-depth gasification) the simulations would develop a large 

deformation ("belly") at the free surface. For carefully selected model parameters and 

procedures, comparisons of the model results and NIST data show favorable agreement. 

However, the research effort was un-able to provide satisfactory explanations for the need to 

significantly increase the model parameters for the in-depth absorption and in-depth gasification 

models or the presence of the "belly" at the free surface of the thermoplastic object. These issues 

remain un-resolved though some efforts continue at NIST to address these issues. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects iii 

The source code, documentation, test problems (all needed files) and presentation material have 

been delivered to NIST. 

 

In the future, it may be possible to couple models such as the condensed phase CFD code to the 

NIST Fire Dynamics Simulator (FDS) code, a CFD model that solves for gas phase transport and 

combustion. The combined tools would enable high-fidelity simulations for fire spread scenarios 

in office environments populated with desktop/deskside personal computers. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects iv 

TABLE OF CONTENTS 

 Page 

DISCLAIMER................................................................................................................................. i 

EXECUTIVE SUMMARY ............................................................................................................ ii

TABLE OF CONTENTS............................................................................................................... iv 

LIST OF FIGURES ....................................................................................................................... vi 

LIST OF TABLES....................................................................................................................... viii 

ACRONYMS AND ABBREVIATIONS...................................................................................... ix 

1. BACKGROUND ........................................................................................................................ 1 

1.1 Introduction................................................................................................................. 1 

1.2 Project Management - Statement of Work, Deliverables and Reporting.................... 2 

1.3 Report Organization.................................................................................................... 3 

2. METHODS AND APPROACH ............................................................................................... 4 

2.1 Overview of Time Varying Grid CFD Model (CPCFD) model................................. 5 

2.2 Details of CPCFD model ............................................................................................ 7 

2.2.1 Tracking The Free Surface ............................................................................ 7 

2.2.2 Wall Boundary Condition............................................................................ 10 

2.2.3 Time Integration Algorithm......................................................................... 12 

3. RESULTS AND DISCUSSION............................................................................................... 14 

3.1 Task 0.......................................................................................................................... 15 

3.2 Task 1.......................................................................................................................... 18 

3.2.1 Task Overview............................................................................................. 18 

3.2.2 Free Surface Smoothing............................................................................... 19 

3.2.3 Tracking Free Surface at the Edges ............................................................. 21 

3.2.4 Task 1B: Modeling Results of Resin PP702................................................ 25 

3.2.5 Task 1C: Moedling Results of Resin PP6523.............................................. 32 

3.3 Task 2.......................................................................................................................... 34 

3.3.1 Task Overview............................................................................................. 34 

3.3.2 Heat Loss and In-depth Absorption Models................................................ 35 

3.3.3 Computational Grid and Case Summary ..................................................... 37 

3.3.4 Results and Discussion ................................................................................ 38 

3.4 Task 3.......................................................................................................................... 41 

3.4.1 Task Overview............................................................................................. 41 

3.4.2 In-depth Gasification Model........................................................................ 42 

3.4.3 Task 3A – CFD Simulation of Resin PP702N............................................. 43 

3.4.4 Task 3B – CFD Simulation of Resin PP6523.............................................. 50 

3.5 Task 4.......................................................................................................................... 51 

3.5.1 Task Overview............................................................................................. 51 

3.5.2 Conjugate Flow-Heat Transfer Problem...................................................... 53 

3.5.3 Computational Grid, Melt Flow Source, Boundary and Initial Conditions. 55 

3.5.4 Implementation of Surface Tension Model ................................................. 57 

3.5.5 Parametric Study of Grid Sensitivity and Time Accuracy .......................... 59 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects v 

3.5.6 Task 4A, 4B and 4C Results and Discussions............................................. 63 

4. CONCLUSIONS....................................................................................................................... 67 

5. LITERATURE REFERENCES................................................................................................ 69 

APPENDIX A – Detailed PowerPoint Summary of Modeling Results .................................. A-1 

Task 0.................................................................................................A-2 

Task 1................................................................................................ A-9 

Task 2................................................................................................A-83 

Task 3..............................................................................................A-153 

Task 4..............................................................................................A-224 

APPENDIX B – Code Documentation.................................................................................... B-1

1.0 Instructions for using the CPCFD code to solve problems of Task 1, 2 and 3.........B-2 

 2.0 Instructions for using the CPCFD code to solve problems of Task 4..................B-8 

3.0 Instructions for adjusting control parameters to avoid divergence.....................B-15 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects vi 

LIST OF FIGURES 

 Page 

Figure 2.1 Sketch of free surface correction procedure.................................................................. 8 

Figure 2.2 Sketch of interior grid adjustment ................................................................................. 9

Figure 2.3 Sketch of the sheme to enforce adiabatic wall boundary condition ............................ 11 

Figure 3.1 Sketch of 2D test problem used in Task 0................................................................... 15 

Figure 3.2 Object height as a function of time.............................................................................. 16 

Figure 3.3 Temperature distribution at a selected time (1000s) ................................................... 16 

Figure 3.4 Velocity vectors colored by velocity magnitude at a selected time (1000s) ............... 17 

Figure 3.5 Geometry and boundary conditions of the 2D test problem for Task 1 ...................... 18 

Figure 3.6 Surface smoothing scheme.......................................................................................... 19

Figure 3.7 Fluctuations in free surface will small blending factors.............................................. 19 

Figure 3.8 Free surface evolution with large blending factors (0.1,0.5,1.0)................................. 20 

Figure 3.9 Schematic for tracking the free surface at the edge..................................................... 21 

Figure 3.10 Free surface evolution using three different surface tracking algorithms ................. 23 

Figure 3.11 Mass loss percentage as function of time for three surface tracking algorithms....... 23 

Figure 3.12 Viscosity-temperature relationships for PP702......................................................... 25 

Figure 3.13 Computational grids used in Task 1A and Task 1B.................................................. 26 

Figure 3.14 Typical residual histories in one time step ................................................................ 27 

Figure 3.15 Mass loss percentage as function of time on four grids (PP702N) ........................... 28 

Figure 3.16 Free surface location and temperature field as a function of time (PP702N) ........... 29 

Figure 3.17 Velocity vectors at Time = 220 seconds (PP702N)................................................... 30 

Figure 3.18 Temperature profiles on the free surface at selected times (PP702N) ...................... 31 

Figure 3.19 Viscosity-temperature relationships for PP6523....................................................... 32 

Figure 3.20 Mass loss percentage as function of time on four grids (PP6523) ............................ 32 

Figure 3.21 Sketch of in-depth absorption and heat loss models ................................................. 35 

Figure 3.22 Incident heat flux distribution in the horizontal direction......................................... 36 

Figure 3.23 Task 2 computational grid at t = 0 sec....................................................................... 37 

Figure 3.24 Mass loss percentage as a function of time for Task 2 cases .................................... 38 

Figure 3.25 Surface temperature and velocity magnitude profiles along the free surface ........... 39 

Figure 3.26 Velocity vectors in Case 2A and Case 2B................................................................. 40 

Figure 3.27 Sketch of in-depth absorption, gasification, and surface heat loss models ............... 42 

Figure 3.28 Sketch of modeled resin in Task 3 ............................................................................ 43 

Figure 3.29 Normalized incident heat flux profile as a function of distance from free surface ... 44 

Figure 3.30 Impacts of changing gasification and absorption models in Task 3A....................... 45 

Figure 3.31 Impacts of external heat flux in Task 3A .................................................................. 46 

Figure 3.32 Comparison of predicted and measured total mass loss rate, fraction of mass loss 

occuring as melt flow and surface temperature for four cases in Task 3A................ 47 

Figure 3.33 Data for fraction of incident radiation versus depth used to define in-depth 

 absorption model A4................................................................................................. 48 

Figure 3.34 Comparison of predicted and measured total mass loss rate, fraction of 

mass loss occurring as melt flow and the surface temperature for two additional 

cases in Task 3A ........................................................................................................ 49 

Figure 3.35 Sketch of two scenarios of melt flow with a catch surface ....................................... 51 

Figure 3.36 Sketch of Task 4 experiment ..................................................................................... 53 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects vii 

Figure 3.37 Conjugated computational model used in Task 4...................................................... 54 

Figure 3.38 Sketch of computational grids adopted in Task 4 modeling ..................................... 55 

Figure 3.39 1-D heat transfer problem for the catch plate............................................................ 56 

Figure 3.40 Diagram of the alogrithm for calculating surface curvature ..................................... 58 

Figure 3.41 Free surface shapes at selected times for different time-step sizes ........................... 60 

Figure 3.42 Mean spreading velocity as a function of the square of the 

maximum time-step size ........................................................................................... 60 

Figure 3.43 Free surface shapes at selected times using different grid. 

No surface tension effects included........................................................................... 61 

Figure 3.44 Free surface shapes at selected times using different grid. 

Surface tension effects included ................................................................................ 61 

Figure 3.45 Mean spreading velocity as a function of the square of the normalized grid size..... 62 

Figure 3.46 Melt flow rate used in Case 4A and Case 4C............................................................ 64 

Figure 3.47 Melt flow rate used in Case 4B ................................................................................. 64 

Figure 3.48 Temperature contours of the melt pool and catch plate ............................................ 65 

Figure 3.49 Mean velocity of melt front as a function of time..................................................... 66 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects viii 

LIST OF TABLES 

 Page 

Table 1.1 Project Tasks and Status ................................................................................................. 2 

Table 3.1 Clock Time and Memory Requirement for Task 1B.................................................... 31 

Table 3.2 Calculated Steady-State Mass Loss Rate...................................................................... 33 

Table 3.3 Task 2 Case Summary .................................................................................................. 37 

Table 3.4 Calculated Steady-State Mass Loss Rate...................................................................... 38 

Table 3.5 In-Depth Gasification Model Constants in Task 3A .................................................... 43 

Table 3.6 Task 3A Case Summary ............................................................................................... 45

Table 3.7 Task 3A Additional Case Summary ............................................................................. 48 

Table 3.8 In-Depth Gasification Model Constants in Task 3B..................................................... 50 

Table 3.9 Task 3B Case Summary ............................................................................................... 50 

Table 3.10 Task 3B Results Summary ......................................................................................... 50 

Table 3.11 Material Properties Used in Task 4 CFD Model ........................................................ 54 

Table 3.12 Parametric Cases for Studying Grid Sensitivity and Time Accuracy......................... 59 

Table 3.13 Task 4 Case Summary ................................................................................................ 63 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects ix 

ACRONYMS AND ABBREVIATIONS 

 

BFC Body Fitted Coordinate 

CFD Computational Fluid Dynamics 

CPCFD Condensed Phase Computational Fluid Dynamics code 

FDS Fire Dynamics Simulator 

NIST National Institute of Standards and Technology 

PP702N Polypropylene resin 702N 

PP6523 Polypropylene resin 6523 

REI Reaction Engineering International 

SIMPLE Semi-Implicit Method for Pressure Linked Equations 

TVD Total Variation Diminishing 

VOF Volume of Fluid 

1D One Dimensional 

2D Two Dimensional 

3D Three Dimensional 

 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 1 

1. BACKGROUND 

1.1 Introduction 

This report describes the work effort by Reaction Engineering International (REI) to develop, 

demonstrate and deliver to the National Institute of Standards and Technology (NIST) a 

condensed phase computational fluid dynamics (CFD) based tool to model the processes of 

melting, flow and gasification of thermoplastic materials exposed to a high heat flux. This work 

effort is in response to NIST solicitation PWS SB1341-05-C-0041. Potential applications of the 

tool include investigating the behavior of polymer materials commonly used in personal 

computers and computer monitors if exposed to an intense heat flux, such as occurs during a fire 

The condensed phase CFD (CPCFD) model delivered to NIST is based on a time dependent 

(time varying) grid CFD method. The time dependent grid CFD model provides a rigorous 

approach, enabling detailed studies of the melting and melt flow processes for thermoplastic 

materials. The model was developed by re-using portions of an existing CFD tool and enhancing 

it with sub-models required to address the specific needs of modeling thermoplastic flows. The 

time integration, spatial discretization and overall solution procedure are based on standard CFD 

methods from the literature. 

CPCFD is written in FORTRAN 90 in an object-oriented form. A 3D, finite volume, multi-block 

body-fitted time dependent (time varying) grid formulation is used to solve the unsteady Navier 

Stokes equations. A multi-grid method is used to accelerate convergence at each time step. Submodels are included to describe the temperature dependent viscosity relationship and in-depth 

gasification and absorption of thermoplastic materials, free surface flows and surface tension. 

NIST data has been used for key material properties as well as to develop and benchmark the 

sub-models for thermoplastic materials. A variety of boundary conditions can be used for the 

velocity field (no-slip, free-slip) and heat transfer to the object (adiabatic, heat loss, specified 

heat flux). 

Model outputs include the time dependent velocity, temperature and position (displacement) at 

points in the thermoplastic body. The model outputs can be imported to standard CFD 

visualization packages. Depending on the sophistication of the visualization software, either 

static images at a selected point in time or animations that highlight the time-varying behavior of 

the thermoplastic object can be produced. The CFD images contained in this report have been 

generated using FIELDVIEW, a CFD visualization package available from Intelligent Light 

(http://www.ilight.com). Additional model outputs include the time history of the mass loss rate 

and heat fluxes that can be used in comparisons to experimental data. 

The accuracy and capabilities of the modeling tool are demonstrated on a series of test cases of 

increasing complexity. The test cases include grid sensitivity studies, adding heat loss boundary 

conditions and in-depth gasification and adsorption. Simulations for two thermoplastic materials 

(PP702N, PP6523) are performed for different heat flux scenarios and test problem 

configurations. 

All calculations presented in this report have been performed with computational resources 

typical of what is found in an engineering office (i.e., a high end personal computer). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 2 

1.2 Project Management - Statement of Work, Deliverables and Reporting 

Highlighted in Table 1.1 are the tasks that were to be performed to execute this project. Included 

in the table is a summary of the status of each Task. The work effort for the project was divided 

into a Base Period, with a subsequent Option Period. As noted in the table, some tasks in the 

Base Period were deleted from the statement of work by NIST. In addition, NIST elected not to 

fund the tasks in the Option Period. 

Table 1.1. Project Tasks and Status. 

Status 

Task 0: Port CPCFD computer code to Fortran 90 Completed 

Task 1: 2D condensed phase melt-drip model with steady imposed heat 

flux. 

• Task 1A = visit NIST for project kickoff meeting/discussion Completed 

• Task 1B = Perform model for Resin 1 Completed 

• Task 1C = Perform model for Resin 2 Completed 

Task 2: Add radiative and convective heat losses and in-depth absorption 

of radiative heat flux 

• Task 2A = perform model for Resin 1 using heat loss 

boundary condition at the heated surface 

Completed 

• Task 2B = perform model for Resin 1 using heat loss 

boundary condition at the heated surface and in-depth 

absorption 

Completed 

• Task 2C = perform model for Resin 2 using heat loss 

boundary condition at the heated surface and in-depth 

absorption 

Completed 

Task 3: Add gasification 

• Task 3A = Perform model for Resin 1 Completed 

• Task 3B = Perform model for Resin 2 Completed 

Task 4: Include melt pool on catch surface 

• Task 4A = perform model for resin 1 with gap 1 Completed 

• Task 4B = perform model for resin 2 with gap 1 Completed 

• Task 4C = perform model for resin 1 with gap2 Completed 

• Task 4D = NIST Site visit – Project update Deleted from SOW 

by NIST 

• Task 4E = NIST Site visit – Base Period Program Review Deleted from SOW 

by NIST 

• Task 4F = Base Period Summary Report Completed 

Option Period 

Task 5: Provide expertise during coupling of condensed phase model 

with gas phase model 

Option Period Not 

Exercised by NIST 

Task 6: 3D flame spread model Option Period Not 

Exercised by NIST 

Summary Report – Final Option Period Not 

Exercised by NIST 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 3 

Project Kick-off meeting

A project kick-off meeting was held at the start of the project, on November 2, 2005 at the NIST 

offices in Gaithersburg, Maryland, to discuss project plan, technical approach and deliverables. 

Deliverables and Reporting

All contractually required reports, forms, technical updates and source code, documentation and 

supporting files for CPCFD have been provided to the NIST COTR. 

The NIST COTR has been provided with the materials listed below. 

1. A one page monthly progress report for each month of the project. 

2. Detailed PowerPoint summaries containing quantitative results as the work was 

performed (technical updates) and at Task completion. The summaries included 

descriptions of model enhancements and model results (e.g., CFD plots for velocity, 

temperature field, animations of solution variables, computational run time, XY plots for 

key results such as mass loss rate, comparisons to NIST data). 

3. The CPCFD computational model at the completion of each Task; provided were all 

needed source code, input files and output files. 

4. A Final Report (i.e., this document) that includes a description of the modeling method, 

test problems for each Task, quantitative results for each task and difficulties encountered 

during model development. 

During the course of the project, conference call based meetings were conducted that included 

NIST and REI technical personnel in which technical updates and technical details of the 

modeling effort were discussed. The conference call meetings were held on an as-needed basis 

(e.g., one to two conference calls per month). 

1.3 Report Organization 

The remainder of this document is organized as follows: 

• Chapter 2 - a description of the methods and approach used in CPCFD; 

• Chapter 3 – a detailed description of the modeling results for each Task performed during 

the Base Period; 

• Chapter 4 – conclusions from the project; 

• Appendix A – a detailed PowerPoint summary for the work effort, model enhancements 

and model results for this project; portions of the material in Appendix A were provided 

to and discussed with Dr. Kathy Butler (NIST COTR) and Dr. Tom Ohlemiller (NIST) 

during the course of the project; 

• Appendix B – CPCFD code documentation, including information on: 

o Instructions for Using CPCFD to Solve Problems of Task 1, 2 and 3; 

o Instructions for Using CPCFD to Solve Problems of Task 4; and 

o Instructions for Adjusting Control Parameters to Avoid Divergence. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 4 

2. METHODS AND APPROACH 

The general approach used in CPCFD is based on standard finite volume CFD methods, the 

details of which are available in many references and thus only briefly described here. Two 

topics that merit a more detailed discussion are the technique used to generate and smooth the 

free surface of the time varying grid and the algorithm used to advance the solution through time 

in a time-accurate manner. 

In the following sections, a description of the numerical details used within CPCFD are 

described. Provided are: 

• an overview of the time varying grid CFD model; 

• details for: 

o tracking the free surface; 

o wall boundary conditions; and 

o the time integration algorithm. 

The description of the methods assumes the reader is familiar with CFD modeling techniques. 

Literature references are noted where appropriate. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 5 

2.1 Overview of Time Varying Grid CFD Model (CPCFD) model 

A numerical method featuring a time-varying grid system was employed to solve the targeted 

thermoplastic melting problems. In the method, a grid with a body fitted coordinate (BFC) 

system was used to discretize the thermoplastic object. As the object deforms the grid deforms 

accordingly, i.e. the free surface of the object is always a boundary of the computational domain. 

The coordinate transformation for a three-dimensional case can be written as 

t

x x x t

x x x t

x x x t

=

=

=

=

τ

ξ ξ

ξ ξ

ξ ξ

( , , , )

( , , , )

( , , , )

3 3 1 2 3

2 2 1 2 3

1 1 1 2 3

 (1) 

where t is time, ) ( , , 1 2 3 x x x is another notation for the Cartesian coordinates (x, y,z), and τ is 

the time in the curvilinear coordinate system ) ( , , ξ 1 ξ 2 ξ 3 . Note thatξ 1 , ξ 2 and ξ 3 are functions of 

time. The Navier-Stokes equations in the new coordinate system is written as 

( )

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂

∂

∂

∂

∂ =⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂

+

∂

∂

∂

∂

+

∂

∂

+⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂ ⎥ =

⎦

⎤ ⎢

⎣

⎡

∂

∂

+

∂

∂

+

∂

∂

∂

∂

+

∂

∂

=⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

∂

∂

+

∂

∂

∂

∂

+

∂

∂

k

j

k

m

k j m

j

p k

j

k

k

p p

i

k

j

ik

k j

j

i k ij

j

k

k

i i

k

j

k

j

k

k

x J x

T k

J x

c Tu

t

c T

J

c T

J

g J x J xu u p

t

u

J

u

J

J x

u

J J t

ξ ξ

ξ ξ

ξ

ρ ξ

ξ

ξ

ρ

τ

ρ

ρ

ξ

σ

ξ

ξ

ρ δ

ξ

ξ

ξ

ρ

τ

ρ

ξ

ρ ξ

ξ

ξ

ρ

τ

ρ

1 1 1 1

1 1 1 1

0 1 1 1

 (2) 

where 

( , , )

( , , )

1 2 3

1 2 3

x x x

J

∂

∂ = ξ ξ ξ is the Jacobian of the transformation, and ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂

+

∂

∂ =

i

k

k

i

ik

x

u

x

u

σ μ is the 

viscous stress tensor [5]. The second term at the left-hand-side of each of three equations above 

comes from the time-dependent coordinate transformation. It represents the additional fluxes of 

mass, momentum and enthalpy (for continuity, momentum and energy equations, respectively) 

introduced by the moving coordinates. It should be noted that the coordinate transformation 

(time varying grid) introduces no approximation to the NS equations. 

A pressure correction method [6] was used to solve the governing equations. All equations are 

discretized using a finite volume method. For the momentum equations the second-order central 

difference scheme is used for spatial discretization due to the large viscosity that occurs in 

thermoplastic flow problems [10]. For the energy equation a second-order Total Variation 

Diminishing (TVD) scheme is applied [10]. Time integration is performed with a fully implicit 

three-level method that gives second-order temporal accuracy [5]. The solution must be tightly 

converged at each time step to maintain numerical stability over long time periods. Hence, a 

multi-grid algorithm is applied to accelerate the convergence of the solution at each time step [5], 

[10]. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 6 

The solution algorithm follows the well-known Semi-Implicit Method for Pressure Linked 

Equations (SIMPLE) approach [6]. In the SIMPLE the solution is advanced through time using a 

series of inner iteration loops and outer iteration loops. For a moving grid where the time history 

of the grid movement is prescribed (e.g., a piston moving in a chamber), the SIMPLE scheme 

can be described as follows: 

• Step 1: All variables are assigned initial values at 0 t = t . 

• Step 2: Time is advanced to t = t + ∆t 1 0 , the grid is moved to the new position and an iterative 

procedure is performed to determine the solution of the coupled non-linear equations at the 

new time level. 

(Note: for the first time step, a first-order Euler backward scheme is used for the time 

integration instead of the three-level method which requires information at two previous time 

levels). 

3⁄4 Within each time step, the following steps are repeated: 

 (Steps 2a to 2d constitute one outer iteration loop) 

2a The momentum equations are discretized and linearized, leading to an algebraic 

equation system for each velocity component. These linear equation systems are 

solved iteratively in turn to obtain an improved estimate of the velocity at the new 

time level. Here, a multi-grid algorithm is used to accelerate the convergence. The 

iterations in the linear equation solver are called inner iterations

2b The improved velocity field is used to calculate new mass fluxes through the 

control volume faces and to invoke the mass-conservation equation; the result is 

the pressure-correction equation. Upon solving for the pressure correction, the 

mass fluxes, cell-center velocities, and pressure are corrected. 

2c The energy transport equation is solved in the same manner as per step (a) to 

obtain better estimates for the temperature. 

2d Fluid properties (e.g., viscosity) and source terms (e.g., gasification mass loss 

rate) are updated. 

3⁄4 The outer iteration loop (steps 2a to 2d) is repeated until the residual level computed before 

the first inner iteration in each equation becomes sufficiently small. 

When the non-linear, coupled equations are satisfied to a desired tolerance, the time is advanced 

by another∆t , and the above described process (Steps 1-2) repeated. The solution from the 

previous time step serves as the initial guess for the solution at the next time level. 

For the thermoplastic melt flow problems of interest in this report, the grid does not move in a 

prescribed manner. Instead, the grid movement if calculated as part of the solution. The arbitrary 

grid movement complicates the solution procedure. Furthermore, as described in this report, 

special grid smoothing procedures must be performed to avoid simulation divergence. The 

SIMPLE solution procedure used for the melt flow problems performed in this report is 

described in the next section. 

 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 7 

2.2 Details of CPCFD model 

2.2.1 Tracking The Free Surface 

The following boundary conditions have to be satisfied at the free surface: 

• kinematic condition: this boundary condition implies that the free surface is an interface 

between two fluids with no flow through it, i.e.: 

m& fs = 0 (3) 

where “fs” denotes the free surface. 

• dynamic condition: this boundary condition implies that the forces acting on the fluid at the 

free surface are in equilibrium. If the surface tension force and the viscous stresses at the 

free surface are neglected, the dynamic condition reduces to the statement that the pressure 

on both sides of the free surface must be the same. 

In CPCFD, the dynamic boundary condition is implemented by treating the free surface as a 

boundary with a prescribed pressure. In the discretized continuity equation the velocity of the 

fluid at the free surface boundary is obtained with a one-sided extrapolation from the interior. 

The pressure-correction equation employs a Dirichlet boundary condition at the free surface. 

With this approach, at the end of a SIMPLE-step the dynamic condition is satisfied - but not the 

kinematic condition. The grid cell face velocities at the free surface, which are corrected after 

solving the pressure-correction equation to enforce local and global mass conservation, can lead 

to non-zero mass fluxes ( mfs & ) To satisfy the kinematic condition, the grid cell vertices which 

define the boundary cell face must be moved so that the correction of the volume swept by the 

free surface (Vfs′ & ) compensates for the mass flux created in the preceding step: 

+ ′ = 0 mfs Vfs

& & (4) 

An iterative procedure for correcting the free surface is used in CPCFD. The procedure is 

incorporated into the outer iteration loop of the SIMPLE algorithm described in Section 2.1. 

The SIMPLE solution algorithm, including free surface tracking, used to perform the melt flow 

simulations in this report is summarized below. 

• Step 1: All variables are assigned initial values at 0 t = t . 

• Step 2: Time is advanced to t = t + ∆t 1 0 , the grid is moved to the new position and an iterative 

procedure is performed to determine the solution of the coupled non-linear equations at the 

new time level. 

(Note: for the first time step, a first-order Euler backward scheme is used for the time 

integration instead of the three-level method which requires information at two previous time 

levels). 

3⁄4 Within each time step, the following steps are repeated: 

 (Steps 2a to 2e constitute one outer iteration loop) 

2a Solve the momentum equations using the geometry defined by the current shape 

of the free surface and the prescribed pressure on it. This replaces the inner 

iteration loop (step 2a) outlined in Section 2.1. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 8 

2b Enforce local mass conservation in each control volume by solving the pressurecorrection equation, using the prescribed pressure boundary condition at the 

current free surface. Mass is conserved both globally and in each control volume, 

but non-zero mass fluxes through the free surface can result. 

2c Correct the position of the free surface so that the volume defined by its corrected 

and previous position compensates for the mass fluxes through the free surface 

obtained in the preceding step. 

2d The energy transport equation is solved in the same manner as per step 2a to 

obtain better estimates for the temperature. 

2e Fluid properties (e.g., viscosity) and source terms (e.g., gasification mass loss 

rate) are updated. 

3⁄4 The outer iteration loop (steps 2a-2e) are repeated until all equations and boundary 

conditions are satisfied (i.e., until the mass fluxes are zero and the residual of each 

equation is smaller than a threshold value). 

The time is advanced by another∆t , and the above described process (Steps 1-2) repeated. The 

solution from the previous time step serves as the initial guess for the solution at the next time 

level. 

The free surface correction procedure described in step 2c (above) is illustrated in Figure 2.1. 

In the current implementation, the vertices that define the free surface are only allowed to move 

in one direction (i.e., in the X direction in Figure 2.1). This significantly simplifies the reVertices at free surface before adjustment 

Control points at free surface before adjustment 

C1 V1 C2

∆hC2

∆hC1

X

Y 

∆hV1

Vertices at free surface after adjustment 

Control points after adjustment

Grid lines before adjustment

Grid lines after adjustment

V0 V2

mfs1 & mfs2 &

Figure 2.1 Sketch for free surface correction procedure. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 9 

meshing algorithm and has been shown to be adequate for modeling the targeted thermoplastic 

deformation problems. 

• The first step in the free surface correction algorithm is to define a set of control points 

along the current free surface. The current cell face centers are selected for this purpose. 

The corrected positions of the control points defined by the displacement ∆hC are 

calculated as 

A

m t

h fs

C

ρ

∆

∆ = & , (5) 

where mfs & is the mass flux thought the cell face, ∆t is the time step, ρ is the density, and 

A is the area of the cell face. 

• The corrected positions of vertices are defined by the displacement∆hV , which is 

computed by averaging neighboring ∆hC ’s. If uniform grid spacing is used in the Y 

direction, such as shown in Figure 2.1, then ∆hV is expressed as 

( ) , , , 1 2

1 ∆hV n = ∆hC n + ∆hC n+ , (6) 

where n is the index of the vertex. 

• The free surface after correction is defined by the corrected vertices shown as dashed line 

in Figure 2.1. 

The grid in the interior of the solution domain must respond to the movement of the vertices in 

the free surface, to preserve its good quality (i.e., avoid grid skewness). A simple and general remeshing technique for the internal grid, which can be successfully applied in all situations and 

for all geometries, is difficult to devise. The current algorithm simply re-adjusts the location of 

the inner vertices to preserve the same relative position between the free surface and the opposite 

boundary. In addition to the adjustment of the interior grid, the mass fluxes through the interior 

cell faces may also need adjustment to maintain mass conservation in the new control volume. 

As shown in Figure 2.2, the mass flux k+1

m& after the grid adjustment is calculated as 

t

V

m m

k k

∆ = − + ρ & & 1 . (7) 

 

Figure 2.2 Sketch of interior grid adjustment. 

Interior vertices after adjustment 

Interior vertices before adjustment

Grid lines before adjustment

Grid lines after adjustment

V 

mk 

mk+1
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 10 

2.2.2 Wall Boundary Condition 

In CPCFD, the wall boundary conditions are implemented by using “ghost” cells. Ghost cells are 

grid cells that form the boundary of the computational domain. The cell center quantities such as 

velocity and temperature of a ghost cell are not updated in the procedure of solving the 

governing partial differential equations (PDEs). A “wall cell” is a ghost cell that forms the solid 

wall surrounding the computational domain. A “fluid cell” is a grid cell in the interior of the 

computational domain whose cell center quantities are the solution of the governing PDEs. 

In this project, wall boundary conditions are characterized as no-slip/slip, no-penetration, and 

adiabatic/constant temperature. 

• In CPCFD, the slip wall boundary condition is implemented by equating the cell-center 

velocity components in a wall cell (i.e., ghost cell) with its neighboring fluid cell. 

• The constant temperature condition is implemented by setting the temperature at the wall 

cell (i.e., ghost cell) cell-center to a given temperature. 

• The no-slip, no-penetration wall boundary is implemented by setting the cell-center 

velocity components in the wall cells (i.e., ghost cells) of the boundary to zero. When 

calculating the diffusive source terms in the momentum equations, the velocity 

components at the wall cell face-center are needed. In CPCFD, the velocity components 

at a grid cell face-center are obtained by interpolating the velocity components at the two 

grid cell-centers that share the same grid cell-face. The interpolation factor is determined 

by the respective distances between the grid cell face-center and the two grid cell-centers. 

For the velocity at the wall cell face-center, the interpolation factor is either 0 or 1 

depending on where the wall cell lies relative to the fluid cell. For example, if the grid 

cell above (north of) the fluid cell is a wall cell (i.e., ghost cell) then the interpolation 

factor is 1. With this approach, the velocity components at a wall cell face-center is 

always set equal to the velocity components at the wall cell cell-center, which is always 

zero for a no-slip, no-penetration wall. 

• The adiabatic wall boundary condition is implemented by setting the neighboring 

coefficient associated with the wall cell to zero. For example, the 1D generic transport 

equation (Eq. (10)) can be written as a discretized algebraic equation for cell P as 

follows: 

( ) AETE + AW TW − AE + AW − S P TP = SU , (8) 

where T is the cell center fluid temperature to be solved; AE, AW, and AP = (AE+AW-SP)

are the coefficients; SU and SP are the source terms; and cell E and cell W are the 

neighboring cells of cell P and lie to the east and west of cell P, respectively. 

For illustration purposes, it is assumed the east side of cell P is an adiabatic wall. Thus, 

cell E is a wall cell (i.e., ghost cell) and its cell center temperature TE must be equal to TP. 

Substituting the condition TE = TP into Eq. (8) results in 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 11 

( ) AW TW − AW − SP TP = SU , (9) 

which implies AP=-(AW-SP). This is equivalent to assigning AE to 0 before AP is 

calculated. In CPCFD, Eq. (9) is used to construct the linear equation system for the 

unknowns and is solved iteratively. 

The above implementation for the adiabatic wall boundary condition implies that the zero 

heat flux condition (i.e., adiabatic condition) at the wall face is satisfied along the line 

that connects the fluid cell cell-center and the wall cell (i.e., ghost cell) cell-center. 

However, if the line is not perpendicular to the wall surface, an error is introduced into 

the solution. In CPCFD, this error is eliminated by adjusting the wall cell (i.e., ghost cell) 

cell-center positions during the re-meshing step so that the line connecting a fluid cell and 

its wall cell (i.e., ghost cell) counterpart is always normal to the wall cell cell-face. This is 

illustrated in Figure 2.3. Here, the wall cell cell-face is defined by vertices a and b, and 

point c is the cross point of the wall cell cell-face and the line connecting the two cell 

centers. As long as point c is within face ab, the scheme remains first-order accurate. If 

point c lies outside ab, then grid refinement in the vertical direction is needed to maintain 

first-order accuracy. 

Adiabatic wall surface

Ghost wall cell 

a c b

Figure 2.3 Sketch of the scheme to enforce adiabatic wall boundary condition. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 12 

2.2.3 Time Integration Algorithm 

As noted above, two time integration methods are used in CPCFD. The first order, implicit Euler 

method is used for the first time step of the solution. A second order implicit method is used for 

all subsequent time steps. 

First Order Euler Method

Originally, the first order implicit Euler method was the only time integration method available 

in CPCFD [5]. For illustration purposes, we consider a one dimensional (1D) generic transport 

equation with constant velocity, constant fluid properties, and no source terms: 

2

2

x x

u

t ∂

Γ ∂

+

∂

∂ = − ∂

∂ φ

ρ

φ φ . (10) 

If a central differencing scheme is used for the spatial derivatives, the discretized algebraic 

equations in the implicit Euler method are written as: 

,

1

1

1

1

1

P

n

W i

n

E i

n AP i + A + A = Q +

−

+

+

+ φ φ φ (11) 

where 

( )

( )

( )

.

2

2

;

; 2

; 2

n

P i

P E W

W

E

t

Q

t

A A A

x x

u A

x

T

x

u A

φ ρ

ρ

ρ

ρ

∆ =

∆ = − + +

∆

Γ − ∆ = −

∆ − ∆ =

 (12) 

The above method is unconditionally stable but has a first order truncation error. 

Second Order Implicit Method

To achieve higher time accuracy, a fully implicit scheme of second order accuracy [5] is used in 

CPCFD. By using a quadratic backward approximation in time and central differencing in space, 

the resulting algebraic equation for the 1D generic transport equation is written as: 

1 1

1

1

1

1

2

+ 2 −

−

+

+

+

∆ − ∆

+ + = n

i

n

i

n

W i

n

E i

n

P i t t

A A A φ ρ φ ρ φ φ φ (13) 

where the coefficients AE and AW are the same as in case of the implicit Euler scheme. However, 

the central coefficient now has a stronger contribution from the time derivative: 

( ) t

AP AE AW ∆ = − + +

2

3ρ , (14) 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 13 

and the source term contains contribution from time level tn-1. The implementation of the second 

order scheme on top of the implicit Euler scheme is straightforward. 

Stability and Time Step Size

Although both methods are unconditionally stable, it was found that the solution may diverge if 

the time step is too large. This can occur due to the free surface tracking scheme being an 

explicit procedure (see Eq. (5)). In the current model, the time step is determined by 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ∆ ∆ = ∆

u

x

t min t ,c max (15) 

where ∆tmax is a user supplied maximum time step and c is the Courant number. In simulating the 

targeted thermoplastic deformation problem, if oscillations occur or the solution diverges, the 

user must either reduce the maximum time step or reduce the Courant number. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 14 

3. RESULTS AND DISCUSSION 

In this section the modeling results for the different test cases investigated are provided. The 

material in this section is organized by Task (see Table 1.1). Provided, in order, are descriptions 

of the work effort and example results for: 

• Task 0: Port CPCFD computer code to Fortran 90 

• Task 1: 2D condensed phase melt-drip model with steady imposed heat flux. 

o Task 1A = visit NIST for project kickoff meeting/discussion 

o Task 1B = Perform model for Resin 1 

o Task 1C = Perform model for Resin 2 

• Task 2: Add radiative and convective heat losses and in-depth absorption of radiative 

heat flux 

o Task 2A = perform model for Resin 1 using heat loss boundary condition at the 

heated surface 

o Task 2B = perform model for Resin 1 using heat loss boundary condition at the 

heated surface and in-depth absorption 

o Task 2C = perform model for Resin 2 using heat loss boundary condition at the 

heated surface and in-depth absorption 

• Task 3: Add gasification 

o Task 3A = Perform model for Resin 1 

o Task 3B = Perform model for Resin 2 

• Task 4: Include melt pool on catch surface 

o Task 4A = perform model for resin 1 with gap 1 

o Task 4B = perform model for resin 2 with gap 1 

o Task 4C = perform model for resin 1 with gap 2 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 15 

3.1. Task 0: Port CPCFD to Fortran 90 

The CPCFD code was originally written in C++. To reduce concerns about future enhancements 

and support, NIST required that the delivered computational model be written in Fortran 90. 

Hence, the first task performed by REI was to re-write the code in Fortran 90. The Fortran 90 

version of CPCFD retains the object-oriented design of the C++ code and is organized into two 

subroutines and a main program. The test problem from the NIST RFP is used to demonstrate the 

consistency between the C++ code and the Fortran code. The results are described below. 

The test problem is a triangular object resting in a flat container. The object is 5 cm high and 

makes an angle of 45 o with the container surface. Figure 3.1 below illustrates the right half of 

the problem. A constant heat flux is applied to the surface of the object, and heat losses due to 

radiation and convection are taken into account. Free-slip boundary conditions and adiabatic 

boundary conditions are applied to the plane of symmetry along the left face of the object and 

also to the base and side wall of the container. 

A 102×52 grid is used to model the test problem. A constant time step of 2 seconds is used. The 

convergence criterion for each time step in the simulation is that at least five significant digits 

remain unchanged from iteration to iteration. The simulation was run for 1300 time steps (2600 

seconds in integration time). 

Shown in Figure 3.2 is the object height as a function of time. The result from the C++ code is 

plotted as a solid line and the result from the Fortran code is plotted as a dotted line. The 

difference between the two solutions are negligible. It requires 2502 seconds for the object 

height, H, to be reduced down to 26Hs 0 (0.013m), after which the ob . ject height decreases at a 

smaller rate as it approaches steady state. Shown in Figure 3.3 and Figure 3.4 are the temperature 

contours and velocity vectors, respectively, at 1000 s. Again, the results are almost identical. It 

should be pointed out that the Fortran code runs slightly slower then the C++ code. 

Figure 3.1 Sketch of 2D test problem used in Task 0. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 16 

Figure 3.2 Object height as a function of time. 

Figure 3.3 Temperature distribution at selected time (1000 s). 

0

0.01

0.02

0.03

0.04

0.05

0.06

0 500 1000 1500 2000

Time (s)

H (x=0) (m)

C++

Fortran90

Temperature (K) time = 1000 s 

C++ Code 

Fortran90 Code 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 17 

Figure 3.4 Velocity vectors colored by velocity magnitude at selected time (1000 s). 

|V| (m/s) time = 1000 s 

C++ Code 

Fortran90 Code 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 18 

3.2. Task 1: 2D model with steady imposed heat flux 

3.2.1 Task Overview 

In this section, the CPCFD code is used to solve for the flow and temperature field that results 

for an upright 2D thermoplastic slab exposed to a steady heat flux on one side. The Task 1 test 

problem is illustrated Figure 3.5. The slab rests on a fixed boundary surface on the right and 

bottom faces with adiabatic heat properties. The geometry of the catch surface is ignored in this 

task. The viscosity of the slab is strongly dependent on temperature, and the material flows out 

over the bottom boundary when its temperature is sufficiently high. The steady imposed heat 

flux is applied directly to the material surface regardless of its shape. 

The behavior of two polypropylene resins, PP702N and PP6523, are examined. As described 

below these two resins have very different viscosity-temperature relationships. This task was 

performed in three steps: 

• Task 1A = visit NIST for project kickoff meeting/discussion; 

The project kick-off meeting was held at the start of the project, on November 2, 

2005 at the NIST offices in Gaithersburg, Maryland, to discuss project plan, 

technical approach and deliverables. 

• Task 1B = Perform model for Resin PP702; and 

• Task 1C = Perform model for Resin PP6523. 

Provided, in order, in the remainder of this section are: 

• the free surface smoothing scheme used to ensure stability during simulations; 

• improvements to the free surface tracking algorithm used to determine the locations of 

vertices that are common to the free surface and a no-slip wall; 

• model results for Task 1B (for PP702N); and 

• model results for Task 1C. (for PP6523). 

Figure 3.5 Geometry and boundary conditions of the 2D test problem for Task 1. 

Left face of sample: 

Constant heat flux q0 applied to 

surface of flowing material: 

k ∂T/∂s = - q0, where s is in the 

direction normal to the surface 

 

Right face of sample: 

u = 0 (no penetration) 

v = 0 (no-slip) 

∂T/∂x = 0 (adiabatic) 

 

Top, bottom faces of sample 

v = 0 (no penetration) 

u = 0 (no-slip) 

∂T/∂y = 0 (adiabatic) 

Ws = 25 mm; hs = 250 mm 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 19 

3.2.2 Free Surface Smoothing 

Modeling experience highlighted the need to develop a 

grid smoothing scheme that could be applied to the free 

surface to avoid fluctuations in the solution when the free 

surface is almost parallel to the gravitational force. As 

shown in Figure 3.6, the new location of vertex C in the X 

direction after one time step is calculated as follows: 

1. Calculate the new face center locations 

t

A

m

x x

CE

old D

D

new

D = + ∆

ρ

& (16) 

2. Calculate the temporary vertex point locations 

( ) new

D

new

C B x = x + x

2

* 1 (17) 

3. Smooth the moving boundary 

( ) ( ) * * 2

4

1 1 C E

old

A

old

C

new

Cx = −α ⋅ x +α ⋅ x + x + x , (18) 

where α is a blending factor, mD & is the mass flux through the cell face, ∆t is the time step, ρ is 

the density, and ACE is the area of the cell face. Note that in Step 2 and Step 3 the positions of the 

boundary vertex point (e.g., point A in Figure 3.6) are not updated and old

A x is used in Eq. (18). 

The algorithms used to calculate the positions of the boundary vertex are described in the next 

section. 

Figure 3.7 shows the evolution of the free surface in two test cases similar to the Task 1 test 

problem. The value of the blending factor α is set to 0.0 and 0.01, respectively. Furthermore, 

note that for α = 0.0 no smoothing is performed at the free surface. 

 

A 

B 

C 

D 

F 

No-slip adiabatic 

Moving 

Ghost 

E 

x 

y 

Figure 3.6 Surface smoothing 

scheme. 

t = 100 

a=0.0 

a=0.01 

t = 200 t = 210 t = 220 t = 230 

Figure 3.7 Fluctuations in free surface with small blending factors. The top row uses a 

blending factor of 0.0 and the bottom row uses a blending factor of 0.01. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 20 

From Figure 3.7, it can be seen that with no smoothing, or a very small blending factor, 

fluctuations develop along the free surface at a very early stage of the simulation. 

Figure 3.8 shows the results for using three larger blending factors. It is found that, when α > 0.1, 

the free surface becomes very smooth and no fluctuations occur during the calculations. For most 

simulations, a blending factor of 1.0 is used. The blending factor is set within the source code for 

CPCFD. 

t = 300 s t = 400 s t = 500 s t = 600 s 

a=0.1 

a=0.5 a=1.

a=1.

a=0.1 

a=0.5

a=0.1,0.5

a=1.

a=0.1,0.5

a=1.

Figure 3.8 Free surface evolution with large blending factors (0.1, 0.5, 1.0). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 21 

Figure 3.9 Schematic for tracking the 

free surface at the edge. 

A 

B 

C 

D 

F 

No-slip adiabatic wall

Moving boundary

Ghost cell 

E 

x 

y 

3.2.3 Tracking the Free Surface at the Edges 

As outlined in Section 2.2.1, the kinematic condition is used to move the control points (i.e., free 

surface cell face centers). The cell vertices between these control points are then moved to 

locations obtained by interpolating coordinates of neighbor control points, but it remains open 

how to move the points at the edges of the free surface. For points on the free surface the vertex 

moves the same distance as the neighboring control point. However, the edges between free 

surface and the solid wall require special attention. 

For the vertex point on the free surface that resides 

on the solid wall the grid vertex point is allowed to 

move along the wall, but the movement is 

performed in a manner that enforces the “no-slip” 

condition. Simply fixing the location of this grid 

point is not adequate; such an approach results in 

significant grid skewness that leads to divergence of 

the simulation. 

Three models (algorithms) have been investigated 

to move the vertex points at the solid wall-free 

surface vertices. In Figure 3.9, point A is a vertex 

shared by the no-slip wall and the free surface, 

vertices A, C, E, and F lie on the free surface, and 

points B and D are cell face centers that serve as 

control points. 

As outlined in Section 2.2.1, the first step in 

determining the moving boundary after each time 

step is to calculate the new control point locations, 

for example 

t

A

m

x x

AC

old B

B

new

B = + ∆

ρ

& (19) 

where mB & is the mass flux through face AC obtained by solving the momentum equations and 

pressure-correction equation. It should be noted that, in solving for mB & , the no-slip boundary 

condition is satisfied. The second step is to calculate new locations of the vertices by 

interpolating the positions of neighbor control points. However, for vertex point A on the edge, 

special treatment is needed. The three models tested pertain to step 2 in the free surface 

correction algorithm in section 2.2.1. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 22 

Model 1

In this algorithm vertex point A is simply moved by the same amount as its neighboring vertex 

(i.e., point C), resulting in: 

( ) new

D

new

C B

new

A x = x = x + x

2

* 1 . (20) 

After computing new

A x , a smoothing step is applied which updates the locations of the interior 

vertices (e.g., point A and point E) and new

A x from Eq. (20) is used instead of oldA x in Eq. (18) . 

No smoothing is applied to new

A x . 

Model 2

The second approach investigated assumes that vertex A moves the same amount as the control 

point next to it (i.e., cell face-center point B). In this approach, Eq. (20) is replaced by 

t

A

m

x x x

AC

old B

B

new

B

new

A = = + ∆

ρ

& . (21) 

After computing new

A x , a smoothing step is applied which updates the locations of the interior 

vertices (e.g., point A and point E) and new

A x from Eq. (21) is used instead of oldA x in Eq. (18) . 

No smoothing is applied to new

A x . 

Model 3

The last approach investigated involves applying the smoothing process to the interior vertices 

before updating the position of the edge vertex. In this approach, Eq. (17) and Eq. (18) are 

applied as described in Section 2.2.1, after which new

A x is computed from: 

( ) * * 2

4

1

C E

old

A

new

C

new

A x = x = x + x + x . (22) 

No smoothing is applied to new

A x . 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 23 

Comparison of the Three Models

Figure 3.10 shows the test results using these three models. The test case is similar to the Task 1 

problem, where the top boundary is a no-slip wall and the right side boundary is the free surface 

exposed to a constant heat flux. 

 

In Figure 3.10, it can be seen that the edge vertex at the top-right corner retreats much faster in 

Model 1 than in Model 2. This is because in Model 2, the movement of point A is approximated 

by the movement of point B rather than point C (see Figure 3.9). The no-slip wall boundary 

condition has a stronger influence on point B than point C because point B is halfway closer to 

the wall than point C. The different melting behaviors observed in Figure 3.10 for Model 1 and 

Model 3 are a result of the surface smoothing process. 

Figure 3.11 shows the mass loss histories 

obtained from the test case using the three 

different surface tracking models 

described above. Although the different 

models strongly impact the shape of the 

free surface (i.e., Figure 3.10), the 

differences in the mass loss histories are 

very small for these simulations. 

Model 1 

Model 3 

Model 2 

t = 250 sec 

Model 1 

Model 2

Model 3 

t = 350 sec t = 450 sec

Model 1

Model 2

Model 3

t = 550 sec

Model 1 Model 2

Model 3 

Figure 3.10 Free surface evolution using three different surface tracking algorithms. 

0

10

20

30

40

50

60

70

0 200 400 600 800 1000

Time (s)

Mass Loss Percentage (%)

Model 1

Model 2

Model 3

Figure 3.11 Mass loss percentage as function 

of time for three surface tracking algorithms. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 24 

Model 1 is used for the simulations performed for the Task 1, mainly because the issue of 

tracking the edge vertex was not raised until the simulations for Task 2 were being performed. 

Furthermore, Model 1 was adequate for the problem addressed in Task 1. During the work effort 

for Task 2 it became necessary to develop a better method to treat the vertex point at the edge of 

the free surface. In the parametric study described in this section, the simulation using Model 3 

predicted a free surface shape that, qualitatively, better agreed with observed behavior in NIST 

experiments and with simulation results obtained by NIST using a volume of fluid (VOF) 

method [1]. Therefore Model 3 is used for the simulations presented for Task 2 and Task 3. None 

of the models were needed for the simulations performed for Task 4. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 25 

3.2.4 Task 1B: Modeling Results for Resin PP702 

 

Figure 3.12 shows the measured viscosity of the 

PP702 polymer as a function of temperature with a 

fixed shear rate of 0.1 [13]. Note that in the figure the 

melt viscosity uses a log scale. In addition, the figure 

highlights a highly non-linear temperature-viscosity 

relationship; that is, small changes in temperature 

result in very large changes in viscosity. The solid line 

in Figure 3.12 is used in the CFD model to describe 

the highly non-linear viscosity-temperature 

relationship, which can be expressed as 

The width and height of the 2D object are 25 mm and 250 mm, respectively. All material 

properties are constant except the viscosity; the density (ρ) is 900 kg/m3

, the heat capacity (Cp) is 

2400 J/kg-K, and the heat conductivity (k ) is 0.5W/m-K. The initial temperature is 300 K and 

the external heat flux is 20,000 W/m2

. The external heat flux is implemented as a source term in 

the energy equation, and it only affects the control volumes immediately next to the free surface 

(i.e., in-depth absorption is not modeled and all the energy is absorbed in the free surface). No 

heat loss is considered in this task. Due to the lack of key physical models such as radiative heat 

loss in the model, overheating of the material (i.e., the temperature rises to more than 1000 K) 

occurs when the width of the object becomes too small). Overheating causes sharp changes in 

viscosity (e.g., one order of magnitude change in viscosity due to several degrees difference in 

temperature) along the free surface. The material with low viscosity tends to roll-over the 

material with high viscosity, which cannot be handled by CPCFD. To overcome this problem, an 

artificial heat sink term is added to the energy transport equation. The heat sink term is 

constructed to mimic the gasification heat loss term. However, no mass loss is involved in the 

model for this Task. The heat sink term can be expressed as a function of temperature (T ) as : 

⎟

⎠

⎞ ⎜

⎝

⎛ = − − T

E

q B exp . (23) 

For Task 1B and Task 1C, the constants in Eq. (23) are chosen as 20 B = 6.276×10 and 

E = −25350. The heat sink term is negligible when the temperature is below 700 K and becomes 

significant when the temperature exceeds 750 K. 

To evaluate the sensitivity of the model results (e.g., mass loss rate) to grid resolution, four grids 

- 14 × 200 , 15× 200 , 19 × 200 and 21× 200 - with different levels of grid resolution near the free 

surface are used in the calculations. The latter three grids are shown in Figure 3.13. Note that the 

lines shown in the figures connect cell centers instead of vertices. In the vertical direction the 

Figure 3.12 Viscosity-temperature 

relationships for PP702. 

( )( ) ( )

( )

( )

( ) T C

C T C

C T C

C T C

T C

f

f T

f T

T f

≥ °

° ≤ < °

° ≤ < °

° ≤ < °

< °

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

− − +

=

425

350 425

200 350

25 200

25

425

10 200 / 200 25 200

10

2

2

1

1

6

6

μ

( ) ( ) 4 2 7 3 f 1 T 10 ^ 14 .48 0.13858 T 5.5960 10 T 7.8665 10 T − − = − + × − ×

( ) ( ) 4 2 f 2 T 10 ^ 53.19 0.2542 T 2.9879 10 T − = − + ×
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 26 

grid cell spacing is uniform with 200 cells. In the horizontal direction the grid cell spacing is 

biased (i.e., non-uniform). In the 15× 200 grid, the first three layers of cells next to the free 

surface have a cell size of 0.625 mm in the horizontal direction. In the19 × 200 grid and the 

21× 200 grid, the cell size near the surface reduces to 0.5 mm and 0.3125 mm, and the number of 

cells with these cell sizes increases to 4 and 6, respectively. These small cells form a layer along 

the free surface with a width of about 2 mm, which is adequate to cover the melt flow based on 

the experience gained in a trial-and-error process. The initial width of the outlet boundary is set 

to a small value, i.e., 0.5 mm. There are 3, 4, and 6 cells in the outlet boundary in the three grids, 

respectively. The ratio of the cell sizes in the horizontal direction is fixed throughout the 

calculation, thus increasing the number of small cells near the free surface leads to higher 

resolution of the sharp gradient near the surface. Because CPCFD uses a moving grid method, 

the region of high grid resolution tracks the time dependent location of the free surface. 

Because the interface is always a boundary in the interface tracking method implemented in 

CPCFD, high grid resolution is only required near the moving boundary. This is not the case for 

methods such as VOF in which high grid resolution must be used throughout the computational 

domain because the interface lies inside the domain and the location of the interface changes 

with time [1]. Hence, using grids with biased cell spacing with dimensions of 15× 200 , 

19 × 200 , and 21× 200 in the CPCFD simulations provides a grid resolution at the free surface 

that is equivalent to that of grids with a uniform cell spacing with dimensions of, respectively, 

40× 200 , 50× 200, and 80× 200 . Note that VOF simulations typically require a grid with a 

uniform grid cell spacing. 

No-slip wall

• 19X200 cells 

• Biased Grid Spacing in 

Horizontal Direction 

9 smallest grid cell 

less than 0.5 mm 

9 four grid cells at 

outlet boundary 

9 Initial grid 

equivalent to 

50X200 uniform 

grid in VOF 

5.e-4m

• 21X200 cells 

• Biased Grid Spacing in 

Horizontal Direction 

9 smallest grid cell 

less than 0.3125 

mm 

9 six grid cells at 

outlet boundary 

9 Initial grid 

equivalent to 

80X200 uniform 

grid in VOF

• 15X200 cells 

• Biased Grid Spacing in 

Horizontal Direction 

9 smallest grid cell 

less than 0.625 

mm 

9 three grid cells at 

outlet boundary 

9 Initial grid 

equivalent to 

40X200 uniform 

grid in VOF

Width of the outlet boundary 

Figure 3.13 Computational grids used in Task 1A and Task 1B. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 27 

In the calculation, the residual is used to determine convergence. The residual for a solved 

quantity Y is defined by the coefficients and source terms of the discretized equation (i.e., AE, AW, 

AN, AS, and Su) and the solved quantities at the cell centers (i.e., YP) and at each cell face center 

(i.e., YE, YW, YN and YS) : 

⎟

⎟

⎟

⎠

⎞

⎜

⎜

⎜

⎝

⎛

Ψ ≡ + ∑

∑ Ψ

Calculation domain

Y

Calculation domain

Y

Y

R

R

0

10 7 log , (24) 

where Y E E W W N N S S U APYP R ≡ A Y + A Y + A Y + A Y + S − 0 , and ΨY ≡ APYP . During each time step, 

iteration is terminated when residuals of all solved quantities are below 1 (i.e., at least six 

significant digits of the solution remain unchanged from iteration to iteration). Figure 3.14 shows 

a typical residual history for a single time step. The number of iterations required to reach 

convergence depends on many factors. Typically, during the early stages of a simulation many 

iterations (e.g., more than 200) are required to converge the solution at each time step. After the 

resin starts to flow, typically less than ten iterations are required. A minimum number of 

iterations in each time step can be set in the input data file to force CPCFD to always perform a 

minimum number of iterations at each time step regardless of the convergence criteria. Other 

factors that impact the convergence history include the under-relaxation factor, time step size, 

and characteristics of the problem itself (e.g., type of resin, external heat flux). 

-6

-5

-4

-3

-2

-1

0

1

2

3

4

5

0 50 100 150 200 250

Number of Iteration

Residual

Mass

U Velocity

V Velocity

Temperature

Figure 3.14 Typical residual histories in one time step. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 28 

Model Results

The calculated mass loss (i.e., the mass leaving the computational domain) percentage as a 

function of time is shown in Figure 3.15 for each of the four grids described above. 

It can be seen that the results are very sensitive to the grid resolution near the free surface. The 

two coarser grids (14x200, 15x200) are not adequate to resolve the sharp temperature gradient 

near the surface and lead to very different melt behaviors. The two finer grids (19x200, 21x200) 

yield similar results in terms of mass loss history and “steady state” mass loss rate (i.e., grid 

convergence). In CPCFD, the size of the time-step is determined based on the CFL constraint, 

which is a model input. The above four cases all used the same CFL number, resulting in each 

case using a different time step. In general, the 21× 200 grid case has the smallest time step size. 

Hence the temporal error also plays a role in making the mass loss histories different. 

Figure 3.16 shows the calculated free surface location and the temperature field as a function of 

time for the case where a 21× 200 grid is used. Within the first 100 seconds, the movement of the 

free surface is negligible under the simulated condition (i.e., 30 KW/m2 heat flux). However, the 

temperature in the thin layer neighboring the free surface can be seen to build up due to the 

external heat flux. As the temperature rises, the viscosity of the resin decreases rapidly. At about 

120 seconds, noticeable movement of the free surface can be seen at the top-right corner. At 

about 200 seconds, the melt material starts to flow rapidly downwards along the free surface. A 

thin layer of melted polymer with high temperature can be seen near the free surface. From about 

200 to 250 second, the mass loss rate appears to reach steady state. After 260 seconds, the flow 

slows down as most of the material has now left the computational domain. 

Figure 3.15 Mass loss percentage as a function of time on four grids (PP702N). 

0

10

20

30

40

50

60

70

80

90

100

0 50 100 150 200 250 300 350 400

Time (s)

Mass Loss Percentage (%)

PP702, 21X200

PP702, 19X200

PP702, 15X200

PP702, 14X200
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 29 

 

In Figure 3.16, it can be seen that the top of the free surface retreats much faster than the other 

part of the surface once the thin melt layer is formed. Between 200s and 215s, a “belly” (or 

bulge) develops in the lower part of the free surface. These phenomena have not been observed 

in experiments or VOF simulations performed by NIST [13]. The cause for this predicted 

behavior is not known. However, parametric studies indicate that: 

1. the free surface smoothing algorithm and the free surface edge tracking algorithm 

described in the previous two sections (i.e., Section 3.2.2 and Section 3.2.3, respectively) 

have great impacts on the shape of the free surface, particularly near the boundary at the 

top wall; and 

2. the inclusion of more physical processes, such as in-depth absorption (see Section 3.3.2 

for details), can help avoid the development of the “belly” that is observed in the 

simulations. 

 

Snap shots of the velocity field (displayed with velocity vectors) are shown in Figure 3.17, which 

shows the melt flow along the free surface. The results shown are from the 21× 200 case and are 

taken at time 220 seconds. Note that the velocity magnitude is zero everywhere except near the 

free surface. Note that the grid lines in the figure connect the cell centers and the vectors are 

formed by using cell center velocity components. In the figure, all velocity vectors are scaled by 

the same factor and the grid cell size varies along the extent of the free surface due to the moving 

grid algorithm used in the solution procedure. 

T=10s T=50s T=100s T=150s T=175s T=200s T=215s T=230s T=240s T=260s 

 Grid: 21X200 

Temperature 

• Thin layer of 

melted polymer 

with high 

temperature

Figure 3.16 Free surface location and temperature field as function of time (PP702N). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 30 

PP702 

Grid: 21X200

Velocity ~ 0 everywhere 

except near free surface layer

Figure 3.17 Velocity vectors at time = 220 seconds (PP702N). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 31 

Compare Model Results and NIST Experimental Data

The temperature profiles on the free surface from the above case (21x200 grid) at selected times 

(200 seconds, 230 seconds, 240 seconds) are plotted in Figure 3.18. The predicted results show 

that, at the early stage when conduction is the main heat transfer mode inside the material, the 

surface temperature distribution is more uniform than that at a later stage, when convective heat 

transfer becomes stronger. The average free surface temperature observed in similar experiments 

at NIST was about 613 K (Ohlmiller FAX 01-18-06). The predicted average surface 

temperatures as shown in the figure are in the range from 675 K to 700 K. The higher average 

surface temperature calculated by the code may be due to the lack of critical models, such as 

radiative heat loss on the surface, gasification, etc. 

 

Computational Resources

The computational times and memory requirements for the cases performed in Task 1B are 

summarized in Table 3.1. All calculations were run on a Dell Precision 450 Desktop PC with a 

Xeon 3 GHz CPU. 

Table 3.1 Clock Time and Memory Requirement for Task 1B Cases. 

Case Wall Clock Time (hr) Memory Requirement (MB) 

14 X 200 24 1 

15 X 200 36 1 

19 X 200 100 2 

21 X 200 120 2 

Figure 3.18 Temperature profiles on the free surface at selected times (PP702N). 

0

0.05

0.1

0.15

0.2

0.25

600 650 700 750

Temperature (K)

Elevation (m)

200 sec (~10% mass loss)

230 sec (~40% mass loss)

240 sec (~60% mass loss)
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 32 

3.2.5 Task 1C: Modeling Results of Resin PP6523 

Figure 3.19 shows the measured viscosity of PP6523 polymer as a function of temperature with a 

fixed shear rate of 0.1 [13]. The solid line in Figure 3.19 is used in the CFD model to describe 

the highly non-linear viscosity-temperature relationship, which can be expressed as 

Model Results

Comparing the viscosity-temperature relationship used for PP6523 with that for PP702N, it can 

be seen that the viscosity of PP6523 is much greater than that of PP702N when the temperature 

is below 620 K, but is an order of magnitude smaller than that of PP702N when temperature is 

above 700 K. 

The grids used in Task 1b are re-used in the current task. Figure 3.20 shows the calculated mass 

loss histories. The case with the coarsest grid (14x200) predicts a very different melting behavior 

as compared to the other three grids; a much longer time is required before the polymer flows 

and the mass loss rate is also much smaller. The differences among the three cases with finer 

grids are much smaller than the 

differences between them and the 

case with the coarsest grid, 

indicating grid convergence has 

been achieved. Again, because the 

time-step size varies differently in 

each case, temporal error also plays 

a role in making the mass loss 

histories different. 

Figure 3.19 Viscosity-temperature 

relationships for PP6523. 

0

10

20

30

40

50

60

70

80

90

100

0 50 100 150 200 250 300 350 400 450

Time (sec)

Percentage of Mass Loss (%)

PP6523, 21X200

PP6523, 19X200

PP6523, 15X200

PP6523, 14X200

Figure 3.20 Mass loss percentage as a function of 

time for PP6523 on four grids. 

( )( ) ( )

( )

( )

( ) T C

C T C

C T C

C T C

T C

f

f T

f T

T f

≥ °

° ≤ < °

° ≤ < °

° ≤ < °

< °

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

− − +

=

430

350 430

220 350

25 220

25

430

10 220 / 220 25 220

10

2

2

1

1

6

6

μ

( ) ( ) 3 2 5 3 8 4

1 f T 10^ 85.843 1.3870 T 7.979 10 T 2.0286 10 T 1.9344 10 T − − − = − + − × + × − ×

( ) ( ) 4 2

2 f T 10^ 100.90 0.4771T 5.533 10 T − = − + ×
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 33 

Compare Model Results for PP702N (Task 1B) and PP6523 (Task 1C)

A “steady-state” mass loss rate (i.e. averaged mass loss rate) is calculated for certain time 

periods and summarized in Table 3.2. For the cases using the finest grids (19x200, 21x200), the 

steady-state mass loss rates are very close to each other. The results also show that PP65223 

flows slightly faster than does PP702N under the same external heat flux, which may be caused 

by the lower viscosity of PP6523 at temperatures above 650K. 

Table 3.2 Calculated Steady-State Mass Loss Rate. 

Case PP702N PP6523 

Mass loss rate 

(g/s m) 

Time period 

(second) 

Mass loss rate 

(g/s m) 

Time period 

(second) 

14X200 28.1 220 - 330 33.8 250 - 400 

15X200 31.5 210 -320 59.1 200 - 290 

19X200 63.3 180 - 260 70.7 200 - 270 

21X200 62.4 180 - 260 70.5 200 - 270 

Computational Resources

The computational time and memory requirement for the cases performed in Task 1C are 

comparable to those of Task 1B (see Table 3.1). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 34 

3.3. Task 2: Surface Heat Loss and In-depth Absorption of Radiative Heat Flux 

3.3.1 Task Overview 

For Task 2, the CPCFD model established in Task 1 is modified to include the radiative and 

convective heat losses from the heated surface and to include in-depth absorption into the 

polymer of the imposed radiative heat flux. Computational results have been obtained for two 

polymer resins (i.e., PP702N and PP6523). The task was performed in three steps 

• Task 2A: perform model for PP702N with heat loss boundary condition at the heated 

surface 

• Task 2B: perform model for PP702N with heat loss boundary condition at the heated 

surface and in-depth absorption 

• Task 2C: perform model for PP6523 with heat loss boundary condition at the heated 

surface and in-depth absorption 

Several modeling issues were investigated during the performance period of Task 2. These 

include the following: 

• a free surface tracking scheme for the no-slip boundary at the top wall; 

• implementation of an adiabatic wall boundary condition for grid lines that are not 

orthogonal to the solid wall; and 

• the free surface smoothing function 

Although not included in the original statement of work, simulation of a third polymer resin 

(PP23K) was also performed to investigate particular aspects of the model. Results and 

discussions for the investigations of the PP23K resin can be found in Section 3.2 and in the 

PowerPoint progress reports and e-mails with NIST dated between May and June, 2006. 

In this section, provided are 

• a description of the surface heat loss model and the in-depth absorption model; and 

• modeling results from Task 2A, 2B, and 2C. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 35 

3.3.2 Heat Loss and In-depth Absorption Models 

The heat loss model added in Task 2 accounts for the heat transfer processes occurring between 

the heated free surface and the surrounding environment. Two heat transfer processes are 

considered - free convection and heat radiation. 

In Task 1, the external heat flux is applied only on the free surface and the energy transmitted 

into the interior of the material occurs only through heat conduction and heat convection. 

However, in experiments, it has been observed that the heated polymer becomes semitransparent and thus allows the external radiative energy to penetrate the exposed surface [13]. 

The in-depth absorption model developed in this task accounts for this mode of energy 

distribution, in addition to heat conduction and heat convection. 

Figure 3.21 is a sketch of the new models developed in this task. A 2D 5.0 cm×10.0 cm (width x 

height) polymer sample is modeled with the free surface facing right. The other three sides of the 

sample are modeled as no-slip adiabatic wall boundaries. The free surface is exposed to an 

external heat flux 

represented by I0, which 

penetrates the surface and its 

intensity (i.e. I) reduces with 

increasing depth due to the 

in-depth absorption effect of 

the material. Heat losses due 

to heat radiation and free 

convection are represented 

by qrad and qconv, 

respectively. These heat loss 

effects are only applied to 

the free surface. Within the 

CPCFD code, qrad and qconv 

are treated as additional 

source terms in the energy 

transport equation and can 

be expressed as: 

( ) 4

0

4 q&radiation = −Aεσ T − T , (25) 

and 

( ) qrconvection = −Ahconv T − T0 & . (26)

In the above equations, A is the surface area; ε is the emissivity (=1.0); σ is the StefanBoltzmann Constant (=5.67e-8 Wm-2K-4); T0 is the ambient temperature (=298 K); and hconv is 

the convective heat transfer coefficient (=8.0 Wm-2K-1). 

The in-depth absorption effect is modeled using an additional heat source term in the energy 

transport equation. The source term can be written as 

Figure 3.21 Sketch of in-depth absorption and 

heatloss models. 

Free surface 

I0

qrad

I

no-slip wall 

Polymer Sample 

(5.0cmX10.0cm) 

qconv

y 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 36 

∫

Ω

= dΩ

dx

dI

qabsorption & , (27) 

where I is the incident heat flux and Ω is the control volume. The model assumes that the indepth absorption effect only happens in the x-direction (i.e., perpendicular to the free surface). 

The magnitude of the incident heat flux inside the material is a function of the distance from the 

free surface in the x-direction and is based on NIST data [13]: 

( )

( )

⎪

⎪

⎩

⎪

⎪

⎨

⎧

= ≤ <

= ≤ <

− − +

− +

10 , 0.3 x 1.3 (cm)

10 , 0 x 0.3 (cm)

2

2

0.5154 1.473 0.4339

0

5.17 6.964

0

x x

x x

I

I

I

I

, (28) 

where I0 is the external heat flux and x represents the distance from the free surface. 

The relationship between the normalized incident heat flux and distance x is plotted in Figure 

3.22. It can be seen that within 2 mm of the free surface about 80% of the energy has been 

absorbed by the polymer, even though the radiative heat continues more than 10 mm into the 

material. Because Eq. (28) does not extend to the total depth of the object modeled (i.e. 5 cm), 

the remaining incident heat flux beyond 1.3 cm is neglected in the current implementation; the 

incident heat flux at this point is less than 2% of the initial external heat flux. 

Figure 3.22 Incident heat flux distribution in the horizontal direction. 

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x, Distance from Free Surface (cm)

Normalized Incident Heat Flux I/I 0
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 37 

3.3.3 Computational Grid and Case Summary 

As noted above, a 5.0 cm×10.0 cm (width x height) 2D block is modeled in Task 2. The 

computational grid is 30 × 80 cells as shown in Figure 3.23. The initial grid resolution is the same 

as the 19 × 200 grid used in Task 

1, and provides resolution 

equivalent to a 100 × 80 grid if 

uniform grid cell spacing is used 

in the horizontal direction and 

the grid cell size is equal to the 

smallest cell near the surface. 

The size of the outlet boundary 

at the bottom-right corner, 

which contains four outlet cells, 

is set to 0.5 mm initially. 

The material properties and 

external heat flux are the same 

as per Task 1. The three cases 

performed in Task 2 are 

summarized in Table 3.3. All 

cases were run on a 3.0 GHz 

Intel Xeon workstation. 

Table 3.3 Task 2 Case Summary. 

Sub-task Resin Case 

Radiative 

and 

convection 

heat losses 

In-depth 

absorption 

Surface 

smoothing 

blending 

factor * 

Mass loss 

and 

simulated 

time period 

Estimated 

wall clock 

time 

Task 2A PP702N Case 2A Included Not 

included 

1.0 60%, 370 s ~50 hrs 

Task 2B PP702N Case 2B Included Included 1.0 65%, 700 s ~26 hrs 

Task 2C PP6523 Case 2C Included Included 0.1 70%, 700 s ~28 hrs 

* See Section 3.2.2 for details 

Figure 3.23 Task 2 computational grid at t = 0 sec. 

No-slip wall

5.e-4m

5.0 cm Free surface

10.0 cm g
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 38 

3.3.4 Results and Discussion 

Plotted in Figure 3.24 are the mass loss histories for the three cases listed in Table 3.3. 

Comparing the results of Case 2A and Case 2B, it can be seen that the inclusion of the in-depth 

absorption model significantly slows down the melting process. The mass loss histories in Case 

2B and Case 2C are about the same. The steady-state mass loss rate of PP6523 from Case 2C is 

slightly larger than that of PP702N from Case 2B, probably due to the lower viscosity of PP6523 

at high temperature. On the other hand, PP702N in Case 2B starts to flow slightly earlier than 

does PP6523 in Case 2C, which may be due to the fact that at low temperatures the viscosity of 

PP6523 is much higher than that of PP702N. 

Listed in Table 3.4 is the “steady-state” mass loss rates calculated from the modeling results. In 

the table, average mass loss rate reported in columns 2-3 are for the time period shown in column 

4. Comparison of the predicted mass loss rate with the experimental data shows that, even with 

the inclusion of the in-depth absorption model, the predicted melt mass flow rate is still much 

higher than the measured data. The lack of other key physical processes in the model such as indepth gasification may explain the discrepancy. 

Table 3.4 Calculated Steady-State Mass Loss Rate. 

Case Mass loss rate (gs-1m-1) Mass loss rate (gs-1) *

 Time period 

Case 2A 24.0 2.4 200 s – 300 s 

Case 2B 6.5 0.65 200 s – 600 s 

Case 2C 7.5 0.75 240 s – 600 s 

Measurement [13] - 0.15 - 

*

 assumes the depth of the object is 0.1m 

0

10

20

30

40

50

60

70

0 100 200 300 400 500 600 700

Time (s)

Mass Loss Percentage (%)

Case 2A

Case 2B

Case 2C

Figure 3.24 Mass loss percentages as a function of time for Task 2 cases. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 39 

Impact of In-Depth Absorption on Model Results

One of the main goals of Task 2 is to assess the impact of the in-depth absorption model on the 

computational results. The in-depth absorption model enables the external energy to be 

distributed within the object beyond the free surface, without resorting to heat conduction and 

heat convection from the surface thus decreasing the temperature gradient near the free surface. 

The impact of including in-depth absorption on the computational results is two-fold. 

1. It reduces the surface temperature and increases the viscosity in the melting layer 

substantially, which in turn slows down the melting flow. This can be seen in Figure 

3.25, which depicts the temperature and velocity profiles along the free surface. All 

profiles are taken at 35% mass loss. Comparison between the temperature profiles in 

Case 2A and Case 2B shows that the average temperature of the free surface drops about 

75 K due to including in-depth absorption in the model. The cooler temperature changes 

the viscosity by at least one order of magnitude, which leads to a significantly different 

surface flow velocity. The velocity magnitude in Case 2A is six times larger than that in 

Case 2B at the same mass loss percentage. In addition, the average surface temperature 

predicted in Case 2B (i.e., about 640 K) is much closer to that observed in similar 

experiments conducted by NIST, which is 648 K (Ohlemiller FAX 01-18-2006). 

Figure 3.25 Surface temperature and velocity magnitude profiles along the free surface at a 

time corresponding to 35% mass loss. 

0

0.02

0.04

0.06

0.08

0.1

500 525 550 575 600 625 650 675 700 725 750

Surface Temperature (K)

Elevation (m)

Case 2A, 255 sec

Case 2B, 382 sec

Case 2C, 365 sec

35% mass loss 35% mass loss 

0

0.02

0.04

0.06

0.08

0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

Surface Velocity Magnitude (m/s)

Elevation (m)

Case 2A, 255 sec

Case 2B, 382 sec

Case 2C, 365 sec
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 40 

2. The second impact of the in-depth absorption model is that it broadens the thickness of 

the melting layer because it distributes radiative energy further into the thermoplastic 

object (i.e., beyond the free surface). This can be seen in Figure 3.26, which shows the 

velocity vectors near the free surface for Case 2A (left) and Case 2B (right) at points in 

time corresponding to (about) 51% and 19% mass loss, respectively. The velocity vectors 

plotted in Figure 3.26 are scaled about the same. The general location of the melting 

layers can be identified by the red vectors in the two large pictures. The four small (inset) 

pictures show that the velocity magnitude in the melting layer for Case 2B is smaller than 

that in Case 2A because the in-depth absorption model reduces the temperature near the 

free surface, which in turn makes the melt flow more viscous. Although both cases 

started with the same initial grid, the grid cell spacing for Case 2A is different than the 

grid cell spacing for Case 2B due to the free surface evolving at a different rate.

 

Figure 3.26 Velocity vectors in Case 2A and Case 2B. 

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-depth absorption) 

~51% mass loss ~19% mass loss 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 41 

3.4. Task 3: In-depth Gasification 

3.4.1 Task Overview

In Task 3, a model describing the in-depth gasification of the thermoplastic material is 

implemented into CPCFD. An empirical Arrhenius expression is used to implement the 

chemistry into the mass and energy equations. Results have been obtained for two polymer 

resins. The task was performed in two steps 

• Task 3A = Perform model for PP702N 

• Task 3B = Perform model for PP6523 

With guidance from the NIST COTR, a parametric study has been performed to select the best 

available rates and parameters for the expressions for the in-depth gasification. Predicted results 

are compared with experimental data provided by NIST. 

In the following are provided: 

• A description of the in-depth gasification model; 

• Model results for resin PP702N (Task 3A); and 

• Model results for resin PP6523 (Task 3B). 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 42 

3.4.2 In-depth Gasification Model

In CPCFD, the gasification process is represented by the removal of polymer mass and energy 

using sink terms in the mass and energy equations. Figure 3.27 is a sketch of the models used in 

the Task 3 simulations. The in-depth gasification model is represented by qgas and mgas, which 

describe the heat and mass losses due to the 

pyrolysis of the polymer melt material. The 

two sink terms are computed as 

⎟ Ω

⎠

⎞ ⎜

⎝

⎛ = − − ∫

Ω

d

RT

E

q& gasification HvρB exp , (29) 

and 

⎟ Ω

⎠

⎞ ⎜

⎝

⎛ = − − ∫

Ω

d

RT

E

m& gasification ρBexp . (30) 

In the above two equations, Hv is the gasification heat (=1250 J/g), ρ is the density, and B and 

E/R are model constants. 

I0

qrad 

I

no-slip wall 

Polymer Sample 

(5.0cmX10.0cm) 

Free surface 

qconv 

y 

Figure 3.27 Sketch of in-depth absorption, 

gasification, and surface heat loss 

models. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 43 

3.4.3 Task 3A – CFD Simulation of Resin PP702N 

The computational domain adopted in Task 3 is shown in Figure 3.28. As per Task 2, a 

10 cm (H)× 5 cm (D) 2D “slice” of the 3D object is modeled. The width of the object (i.e., W = 10 

cm) is only used to compare to test data. The 2D grid is the same as that used in Task 2. 

Multiple gasification models 

have been investigated in order 

to achieve better agreement with 

measured data. The results 

shown in this section use two of 

them. Both models used the same 

functional form, but different 

numerical parameters. The model 

constants are summarized in 

Table 3.5. The mass loss rate due 

to gasification per unit volume y

is related to temperature T as: 

⎟

⎠

⎞ ⎜

⎝

⎛ = − RT

E

y ρB exp . (31) 

Table 3.5 In-Depth Gasification Model Constants in Task 3A. 

 A (s-1) E/R 

Model G1 2.18e+12 24,400 

Model G6 3.27e+14 24,400 

Model G1 is the base model provided by NIST [13]. Model G6 increases the magnitude of the 

gasification rate by a factor of 150. For information on the other in-depth gasification models 

investigated, see Appendix A. 

Description of In-Depth Absorption Models A1 and A3

Several in-depth absorption models have also been used in this Task to study their impacts on the 

gasification process. The results from three of the models are discussed in this section; results 

obtained using the other models can be found in Appendix A. 

1. The first model - Model A1 – was used in Task 2 and is referred to as the base model. 

The normalized incident heat flux can be written as a function of the distance from the 

free surface: 

( )

( )

⎪

⎪

⎩

⎪

⎪

⎨

⎧

= ≤ <

= ≤ <

− − +

− +

10 , 0.3 x 1.3 (cm)

10 , 0 x 0.3 (cm)

2

2

0.5154 1.473 0.4339

0

5.17 6.964

0

x x

x x

I

I

I

I

 (32) 

2D “slice” for 

CFD model

W = 10 cm

D = 5 cm

H = 10 cm

Radiation

Feed Direction

Melting Front

2D “slice” for 

CFD model

W = 10 cm

D = 5 cm

H = 10 cm

Radiation

Feed Direction

Melting Front

W = 10 cm

D = 5 cm

H = 10 cm

Radiation

Feed Direction

Melting Front

Figure 3.28 Sketch of modeled resin in Task 3. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 44 

2. The second model – Model A3– is different from Model A1 in that the incident heat flux 

dissipates much quicker inside the polymer resin. The relationship between the 

normalized incident heat flux and the distance from the free surface is 

exp( 11.513 ), 0 0.4 (cm)

0

= − x ≤ x <

I

I (33) 

The model constant (i.e., -11.513) used in Model A3 is obtained by requiring 99% absorption at 

x = 0.4 cm. This model appears to give the best simulation results when combined with model 

G6 when comparing to experimental data. 

A third model (Model A4) and its simulation results will be described at the end of the section. 

Figure 3.29 shows a comparison of the two in-depth absorption models. It can be seen that the 

base model (i.e., Model A1) allows the incident heat flux to penetrate much further beyond the 

free surface than does Model A3. 

Model Results and Comparison to Experimental Data for Models (G1, A1) and (G6, A3)

Four cases have been performed in Task 3A and are listed in Table 3.6. Radiative and convective 

heat losses at the free surface are included in all cases using the model established in Task 2. 

Case 3A is the baseline case using base models for gasification and absorption. Model G6 and 

Model A3 in Case 3N are developed to match measured mass loss rates at 30 KWm-2 heat flux. 

Case 3P and 3Q are the same as Case 3N but have different heat flux boundary conditions. All 

calculations were run on a Dell Precision 450 Desktop PC with a Xeon 3 GHz CPU. 

Figure 3.29 Normalized incident heat flux profile as a function of distance from free surface. 

0

0.2

0.4

0.6

0.8

1

1.2

0 0.2 0.4 0.6 0.8 1 1.2

Distance (cm)

I/I0

Model A1

Model A3
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 45 

Table 3.6 Task 3A Case Summary. 

Case Gasification 

model 

Absorption 

model 

External heat 

flux (KWm-2) 

Simulated 

time period 

Wall clock 

time 

Case 3A G1 A1 30 550 sec ~24 hrs 

Case 3N G6 A3 30 800 sec ~36 hrs 

Case 3P G6 A3 46 1000 sec ~42 hrs 

Case 3Q G6 A3 20 1000 sec ~42 hrs 

Figure 3.30 shows the mass loss histories for Case 3A and 3N. The two cases use different 

combinations of gasification and absorption models and the results are vastly different. In the 

plot, the solid lines represent the total mass loss (i.e., material leaving the computational domain 

due to gasification and melting). The dashed lines represent the mass loss occurring as melt flow 

only. The difference between solid and dashed line represents the mass loss due to gasification. 

In Case 3A, the fraction of mass loss occurring as melt flow is more than 99% while the mass 

loss due to gasification is negligible. However, the experimental data suggests that, at 30 KWm-2

heat flux, about 30% mass loss is caused by gasification. 

To match the experimental data, the base models in Case 3A were adjusted to form Model A3 

and Model G6. In Model A3, the penetration length of the incident heat flux is significantly 

shortened to increase the temperature near the free surface. In Model G6, the in-depth 

gasification process is intensified by two orders of magnitude. As noted in Table 3.6, Case 3N 

uses these models. From Figure 3.30, it can be seen that for Case 3N that: 1) the overall mass 

loss rate becomes much smaller than that of Case 3A; and 2) the mass loss due to gasification 

becomes more significant. 

Figure 3.30 Impacts of changing gasification and absorption models in Task 3A. 

0

50

100

150

200

250

300

0 200 400 600 800

Time Exposed (s)

Mass Loss (g)

Case 3N

Case 3A 

(fraction of mass loss occurring 

as melt flow > 99%)

0

50

100

150

200

250

300

0 200 400 600 800

Time Exposed (s)

Mass Loss (g)

Case 3N

0

50

100

150

200

250

300

0 200 400 600 800

Time Exposed (s)

Mass Loss (g)

Case 3N

Case 3A 

(fraction of mass loss occurring 

as melt flow > 99%)

solid line = total mass loss. 

dashed line = mass loss due to 

melt flow. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 46 

The impacts of different external heat fluxes on the simulation results are shown in Figure 3.31. 

Three cases, i.e., Case 3P, 3N and 3Q, are compared with each other. Again, the solid lines in the 

plot represent the total mass loss and the dashed lines represent the mass loss occurring as melt 

flow. It can be seen that the overall mass loss rate increases as external heat flux increases. The 

results also suggest that higher external heat flux leads to higher gasification mass loss 

Illustrated in Figure 3.32 are comparisons of the simulation results (i.e., total mass loss rate, 

fraction of mass loss occurring as melt flow, and average free surface temperature) for Case 3A, 

Case 3N, Case 3P and Case 3Q with NIST experimental data. In the figure, the predicted overall 

mass loss rate and the fraction of mass loss occurring as melt flow are calculated from 200 

seconds to the end time of each simulation. The predicted average surface temperatures are 

calculated at 25% total mass loss. 

For the baseline case (Case 3A), the predicted overall mass loss rate and the fraction of mass loss 

occurring as melt flow are too high as compared to the experimental data. After adjusting the 

gasification and absorption models, the agreement between the data and the calculated results 

from Case 3N is satisfactory. Both Case 3A and Case 3N predict the correct average surface 

temperature. 

It is interesting to see that, although Case 3P and 3Q share the same gasification and absorption 

models with Case 3N which are only “optimized” for the 30 KWm-2 data point, they nevertheless 

predict the correct trend shown in the experimental data set. 

0

50

100

150

200

250

0 200 400 600 800

Time (s)

Mass loss (g)

Case 3P 

Case 3N 

Case 3Q 

Figure 3.31 Impacts of external heat flux in Task 3A. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 47 

0.00

0.15

0.30

0.45

0.60

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Total Mass Loss Rate (g/s)

Measurement

Simulation, Case 3N, 3P & 3Q

Simulation, Case 3A

0.0

0.2

0.4

0.6

0.8

1.0

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Fraction of Mass Loss Occuring as Melt Flow

Measurement

Simulation, Case 3N, 3P & 3Q

Simulation, Case 3A

0

100

200

300

400

500

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Surface Temperature (°C)

Measurement

Simulation, Case 3N, 3P & 3Q

Simulation, Case 3A

Figure 3.32 Comparison of predicted and measured total mass loss rate, fraction of mass loss 

occurring as melt flow and the surface temperature for four cases in Task 3A. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 48 

Description of In-Depth Absorption Model A4

A third in-depth absorption model – Model A4 – was used in for simulations. Model A4 is based 

on data obtained in a recent experiment conducted at NIST [13]. The model is shown in Figure 

3.33. A mathematical description for Model A4 is provided in Eq. (34) and shown in Figure 3.33 

where y represents I/I0 in the model formulation. 

0.6878exp( 2.6103 ) 0.3121exp( 0.386 ), 0 0.6 (cm)

0

= − x + − x ≤ x <

I

I (34) 

Model Results and Comparison to Experimental Data for Model (G6, A4)

Two cases are summarized in Table 3.7. Both simulations were performed assuming an external 

heat flux at 30 KWm-2 . 

Table 3.7 Task 3A Additional Case Summary. 

Case Gasification 

model 

Absorption 

model 

External heat 

flux (KWm-2) 

Simulated 

time period 

Wall clock 

time 

Case 3R G6 A4 30 1000 sec ~20 hrs 

Case 3S G1 A4 30 600 sec ~10 hrs 

Illustrated in Figure 3.34 are comparisons of the simulation results (i.e., total mass loss rate, 

fraction of mass loss occurring as melt flow, and average free surface temperature) for Case 3R 

and Case 3S with NIST experimental data. In the figure, the predicted overall mass loss rate and 

the fraction of mass loss occurring as melt flow are calculated from 200 seconds to the end time 

of each simulation. The predicted average surface temperatures are calculated at 25% total mass 

loss. It can be seen that using the original in-depth gasification model G1 (i.e., Case 3S), the 

predicted total mass loss rate is four times higher than the measured value and the predicted 

fraction of mass loss due to gasification is negligible and much smaller than the measured value. 

However, the predicted surface temperature appears to be in good agreement with the data. By 

Model 

Figure 3.33 Data for fraction of incident radiation versus depth used to 

define in-depth absorption model A4. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 49 

increasing the gasification rate by a factor of 150 in Model G6 (i.e., Case 3R), the predicted total 

mass loss rate is 50% higher than the measurement and the predicted fraction of mass loss 

occurring as melt flow and the surface temperature are in good agreement with the data. 

\ 

The results from Case 3A, Case 3P, Case 3R, and Case 3S presented in this section suggest that, 

to match the experimental data, the kinetic parameter for the in-depth gasification in the base 

model (G1) has to be changed significantly. The reason why an adjustment to the in-depth 

gasification model parameter was necessary remains unsolved. 

0.00

0.15

0.30

0.45

0.60

0.75

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Total Mass Loss Rate (g/s)

Measurement

Simulation, Case 3R

Simulation, Case 3S

0.0

0.2

0.4

0.6

0.8

1.0

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Fraction of Mass Loss Occuring as Melt Flow

Measurement

Simulation, Case 3R

Simulation, Case 3S

250

300

350

400

450

500

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Surface Temperature (°C)

Measurement

Simulation, Case 3R

Simulation, Case 3S

Figure 3.34 Comparison of predicted and measured total mass loss rate, fraction of mass 

loss occurring as melt flow and the surface temperature for two additional cases 

in Task 3A.
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 50 

3.4.4 Task 3B – CFD Simulation of Resin PP6523 

Simulations have also been performed for resin PP6523. The model constants in the gasification 

model used in this task are listed in Table 3.8. Model G7 is the base model provided by NIST 

[13]. The magnitude of the gasification rate is increased by a factor of 150 in Model G8 based on 

the experience of modeling PP702N in Task 3A. 

Table 3.8 In-Depth Gasification Model Constants in Task 3B. 

 B (s-1) E/R 

Model G7 2.23e+13 26,000 

Model G8 3.35e+15 26,000 

 

The absorption models are the same as those used in Task 3A. Two cases were run and are 

summarized in Table 3.9. Case 3C is the baseline case. 

Table 3.9 Task 3B Case Summary. 

Case Gasification 

model 

Absorption 

model 

External heat 

flux (KWm-2 

Simulated time 

period (s) 

Wall clock 

time (hr) 

Case 3C G7 A1 30 500 ~24 

Case 3I G8 A3 30 2000 ~60 

The simulation results are summarized in Table 3.10., which are calculated in the same manner 

as those calculated in Task 3A. For the baseline case (Case 3C), behavior that is very similar to 

that of Case 3A is observed; that is, the melting process is very fast and the mass loss due to 

gasification is negligible. In Case 3I, by shortening the heat flux penetration length and 

strengthening the gasification process, the melting process slows down dramatically and more 

than 70% of the overall mass loss is caused by gasification. The baseline case also predicts a 

higher surface temperature (e.g., 20 K higher at 25% mass loss). No experimental data has been 

provided to REI to compare with the model results. 

Table 3.10 Task 3B Results Summary. 

Case Total mass loss rate (gs-1) Fraction of mass loss 

due to melt flow 

Average surface 

temperature at 25% 

mass loss (°C) 

Case 3C 0.69 0.99 382 

Case 3I 0.065 0.28 362 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 51 

3.5. Task 4: Modeling Melt Pool on Catch Surface 

 

3.5.1 Task Overview 

In Task 4, CPCFD is modified to include the flow of melt on the solid catch surface. Included in 

the model are the catch surface material properties and finite thickness, melt solidification and 

surface tension gradient effects. The model also includes a nonzero tilt angle for the catch 

surface. Results are presented for three cases: the two polymer resins with a single specified tilt 

angle of the catch surface, and one of these polymer resins with a different surface tilt angle. The 

Task was performed through the following steps: 

• Task 4A = perform model for resin 1 (PP702N) with angle 1 

• Task 4B = perform model for resin 2 (PP6523) with angle 1 

• Task 4C = perform model for resin 1 (PP702N) with angle 2 

• Task 4D = site visit to NIST by REI 

• Task 4E = site visit to NIST by REI for Base Period project review 

• Task 4F = base period summary report. 

Note that NIST deleted Task 4D and Task 4E from the statement of work and thus are not 

described here. Task 4F (Summary Report) resulted in this document. 

The original plan was to use one grid block (Block 1) to represent the melting slab (as per Tasks 

1-3), a second (Block 2) to represent the overhanging melt flow, and a third (Block 3) to 

represent the flow on the catch surface, as sketched in Figure 3.35. 

After a feasibility study and discussions with NIST, the following decisions were made regarding 

the Task 4 modeling approach. 

1. The simultaneous processes of melt flowing along the free surface of the polymer object and 

the melt flowing along the catch surface should be decoupled and modeled as two separate 

problems. The melt flow on the free surface of the polymer object is to be solved using the 

Figure 3.35 Sketch of two scenarios of melt flow with a catch surface. In schematic at left 

the melt flow spans the gap. In schematic at right the melt flow drips. 

Block 1: Polymer object

Block 2: Overhanging melt

Block 3: Catch basin

Block 1: Polymer object

Block 2: Overhanging melt

Block 3: Catch basin

Melt drop
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 52 

approach established in Task 1, 2 and 3. In Task 4, the focus is to use CPCFD to model the 

melt flow on the catch surface only. 

2. The overhanging melt is a great challenge to CPCFD and will not be included in the CFD 

calculation because it merely provides a bridge to transport the melt flow from the polymer 

object to the catch surface. 

3. The coupling between the simulations of Block 1 and Block 3 is implemented through 

additional source terms in the governing equations. These source terms are computed during 

the process of solving the Block 1 melt problem (i.e., Task 3). These source terms are time 

dependent and are stored in an output/input data file. When solving the Block 3 melt flow 

problem (i.e., Task 4), these time dependent source terms (i.e., mass source, energy source) 

are read into the model at the beginning of the simulation. As the Task 4 simulation 

progresses, the source terms are tabulated and added to the corresponding governing 

equations to mimic the melt flow from the polymer object landing on the catch basin. 

4. The Task 4 simulations also include the heat transfer process occurring within the catch 

plate. The heat transfer and melt flow of the polymer resin on the catch surface and the heat 

transfer inside the catch plate are modeled as a conjugate heat transfer problem. 

In the following sections are provided: 

• a brief description of the conjugate flow-heat transfer problem; 

• a discussion of the surface tension models used in the model of melt flow on the catch 

basin; 

• a grid sensitivity study; 

• a study of the time accuracy of the simulation; and 

• model results and discussions for Task 4A, Task 4B and Task 4C. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 53 

 

3.5.2 Conjugate Flow-Heat Transfer Problem

A sketch of the experimental configuration to be used by NIST is shown in Figure 3.36. In the 

experiment, the hot melt from the polymer slab will flow down and spread on the ceramic catch 

plate. The catch plate has insulated sides and is heated from below by an aluminum plate that is 

attached to the bottom of the ceramic plate. The temperature of the catch surface will be 

sufficient to allow the hot melt to spread without piling up. The temperature distribution along 

the interface between the ceramic plate and the aluminum plate is assumed to be uniform and 

constant due to the large heat conductivity of the aluminum. 

 

At the start of the experiment, the catch plate is heated by the aluminum plate and the 

temperature distribution inside the catch plate is in steady state. After the hot melt (from the test 

piece located above the catch plate) lands on the catch surface, coupled transient heat transfer 

processes occur within the hot melt and the catch plate due to the temperature difference between 

the two materials and the heat losses that occur at the melt pool surface. 

The heat transfer process inside the ceramic catch plate plays a vital role in the formation of the 

melt pool on the catch surface. To calculate the temperature distribution in the catch plate 

accurately, the heat transfer problem in the catch plate and the hot melt flow problem are solved 

simultaneously. A sketch of the computational model is shown in Figure 3.37. It should be noted 

that only a 2D “slice” in the L-H plane is solved in the Task 4 simulations. The width of the 

domain (i.e., W = 7 cm) is used to facilitate the comparisons with data from future experiments. 

In this model, the computational domain includes both the catch plate and the melt pool on the 

catch surface. The bottom of the catch plate is assumed to be a constant temperature boundary. 

The interface between the melt pool and the catch plate lies within the problem domain and does 

not need any special treatment. The upper boundary of the melt pool is treated as a free surface 

with radiative and convective heat losses. 

Ceramic Catch Plate 

Heated Aluminum Plate 

Hot melt from polymer slab 

seeing 30 KWm-2 heat flux 

Sides insulated 

Catch surface 

Interface with uniform and constant 

temperature distribution 

Figure 3.36 Sketch of Task 4 experiment. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 54 

The material properties of the polymer melt and the catch plate used in the CFD calculations are 

listed in Table 3.11. Because the catch plate must remain solid during the calculation, its 

viscosity is set to a very large value to suppress motion in the catch plate. 

Table 3.11 Material Properties Used in Task 4 CFD Model. 

 Polymer melt Catch plate 

Heat capacity (J kg-1 K-1) 2400 961 

Density (kg m-3) 900 2800 

Heat conductivity (wm-1K-1) 0.25 1.26 

Viscosity (Pa-s) Temperature dependent 1.e6 

Figure 3.37 Conjugated computational model used in Task 4. 

Constant T

Catch Surface

Polymer melt

Radiative and Convective Heat Losses

on the free surface of the melt pool

W = 7 cm

L = 24 cm

Catch Plate

Temp. (°C)

H = 1.27 cm

Constant T

Catch Surface

Polymer melt

Radiative and Convective Heat Losses

on the free surface of the melt pool

W = 7 cm

L = 24 cm

Catch Plate

Constant T

Catch Surface

Polymer melt

Radiative and Convective Heat Losses

on the free surface of the melt pool

W = 7 cm

L = 24 cm

Catch Plate

Temp. (°C)

H = 1.27 cm
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 55 

3.5.3 Computational Grid, Melt Flow Source, Boundary and Initial Conditions 

A sketch of the 2D computational grid is shown in Figure 3.38. The grid lines in the picture 

connect cell centers. The grid consists of two parts. 

• The bottom part ( L × H 2 , 24cm ×1.27cm ) uses a fixed grid to model the catch plate. In 

the H2 direction, there are 10 internal cells and the grid cells are stretched with the finest 

grid spacing near the interface between the catch plate and the melt pool. In the L

direction, uniform grid cell spacing is used. 

• The top part ( L × H1 ) uses a moving grid to model the melt pool. In the H1 direction, 

there are 20 internal cells with uniform spacing. The initial value of H1 is 0.1 mm. The 

grid in the L direction uses the same number of cells and the same grid call spacing as the 

L × H 2 part of the grid. 

The melt flow sources (i.e., mass source and energy source), which approximate the melt flow 

from the polymer slab, are added to the discretized governing equations in selected cells at the 

center of the domain in the L direction next to the free surface. The source terms added to the 

momentum and energy equations can be expressed as m(U U ) & flow − and m(T T ) & flow − , 

respectively, where m& is the mass flow rate (i.e., the mass loss rate computed in Task 3), Uflow ( = 

0 ms-1) and Tflow ( = 650 K) are the velocity and temperature of the hot melt, and U and T are the 

solved variables. The melt flow temperature is estimated based on the simulation results from 

Task 3. 

Figure 3.38 Sketch of computational grids adopted in Task 4 modeling. 

L

H2

H1

Melt flow source terms are added in a certain 

number of cells at the center of the domain next to the free surface

free surface

g

L

H2

H1

Melt flow source terms are added in a certain 

number of cells at the center of the domain next to the free surface

free surface

g
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 56 

The top boundary of the melt pool is a free surface with radiative and convective heat losses. The 

free surface is not exposed to any external radiative heat source. The convective heat loss rate 

qconv & is computed as 

( ) qconv = −Ah Ts − Ta & , (35) 

where A is the surface area, h (= 8 Wm-2K-1) is the convection heat transfer coefficient, Ts is the 

surface temperature and Ta (=298 K) is the ambient temperature. The radiative heat loss rate 

qrad & is calculated as 

( ) 4 4

qrad = −A Ts −Ta & εσ , (36) 

where ε is the emissivity (=1.0); σ is the Stefan-Boltzmann Constant (=5.67e-8 Wm-2K-4). 

The bottom boundary of the catch plate is treated as a no-slip constant temperature boundary. For 

resin PP702N the bottom temperature is set to 675 K and for PP6523 the temperature is set to 

700 K. It has been found that these temperatures have a great impact on the simulation. In 

particular, if the bottom temperature is too low, the melt material will pile up on the surface and 

the simulation will diverge. The values used in Task 4 were obtained via a trial-and-error 

process. The two sides of the domain are modeled as outlet boundaries. Here, a zero-gradient 

condition is used on by assigning the cell-centered quantities of the fluid cells next to the outlet 

boundary to the neighboring ghost cells. 

The initial temperature distribution of the catch plate is obtained by solving a 1-D steady state 

heat transfer problem, the results of which are shown in Figure 3.39. 

The energy balance equation for the above problem can be written as: 

( ) ( ) 4 4 Ts Ta Ts Ta h

dy

dT k = − + εσ − , (37) 

where k ( = 1.26 Wm-1K-1) is the heat conductivity of the catch plate, and ∆y (=1.27 cm) is the 

height of the plate. If the bottom temperature TB is 675 K as in modeling PP702N, the initial 

catch plate temperature T can be expressed as 

T = TB − 6860.63(y + 0.0127), - 0.0127 m ≤ y ≤ 0.0 m (38) 

For a 700 K backside temperature as in modeling PP6523, the initial temperature is 

T = TB − 7570.87(y + 0.0127), - 0.0127 m ≤ y ≤ 0.0 m (39) 

The initial cell center temperature for the melt pool is set to the surface temperature Ts computed 

using Eq. (38) and Eq. (39) at y = 0 m. 

Figure 3.39 1-D heat transfer problem for the catch plate. 

TB

TS

Y 

X TS

Qcon Qra TA

∆y Qcon
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 57 

3.5.4 Implementation of Surface Tension Model 

Two surface tension effects have been considered in the current task. 

The first one is the surface tension force normal to the interface that balances the pressure 

difference across the interface. The resulting equation is known as the Young-Laplace equation 

[10]: 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ∆ = +

1 2

1 1

R R

P γ , (40) 

where ∆P is the pressure difference over the interface; R1 and R2 are the principle radii of the 

curvature of the surface; and γ is the surface tension. In the CFD implementation, the pressure 

difference is modeled as a body force in the momentum equations as: 

n

n

n

n

F P A n A n vol

r

r

r

=

= ∇ •

= −∆ ⋅ ⋅ = − ⋅ ⋅ ⋅

ˆ

ˆ

ˆ ˆ

κ

γ κ

 (41) 

where A is the surface area; κ is the curvature of the free surface; and n

r

is the normal vector of 

the free surface. 

The second effect is the so-called surface tension gradient effect in fire spreading studied by 

Sirignano and Glassman [11]. They found that the variation of the surface tension due to the 

temperature gradient along the surface amounts to a stress by which the surface liquid is pulled 

away from the flame front. Thus, the hot liquid is carried forward in the direction of flame 

propagation and energy is transported by convection. At the surface, the stress resulting from the 

variation of the surface tension equals the viscous stress: 

τ γ

τ

γ μ = ≡ ∂∂

d

d

n

U , (42) 

where μ is viscosity, U is the velocity, n is the normal direction and τ is the tangential direction. 

In the CFD implementation, the additional surface shear stress term due to the surface tension 

gradient is modeled as a body force term in the momentum equations as follows: 

τ

τ

τ

γ τ γ γ τ τ τ

v

v

r r r

=

= ⋅ = +

ˆ

S ˆ i j x y

 (43) 

where τ

v is the tangential vector of the free surface. τx γ and τy γ are the additional source terms in 

the u- and v-momentum equations, respectively. Note that the surface tension models are only 

applied to cells that form the free surface (i.e., these source terms are zero for internal cells). 

In CPCFD, the surface tension is correlated to the temperature as [12]: 

γ = −0.04T(°C) + 29.10 (165 - 220 °C) . (44) 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 58 

The correlation has been suggested by NIST [13] and is valid for Eastman Epolene D-10 atactic 

polypropylene. For temperatures outside the range of 165 to 220°C , extrapolation is used in the 

CFD modeling. 

Another important quantity in the surface tension model described above is the surface curvature. 

Figure 3.40 is a diagram of the algorithm developed to calculate the surface curvature κ. Note 

that the current implementation of the algorithm is only valid for the 2D problems solved in Task 

4. The extension of the algorithm to 3D geometries is straightforward. In the algorithm, the 

surface is approximated by line segments that connect vertices. Referring to Figure 3.40, the 

normal vector n1 is calculated based on the vertices (x1, y1) and (x2, y2) and the normal vector n2

is calculated based on the vertices (x2, y2) and (x3, y3). The curvature at vertex (x2, y2) can then be 

approximated as: 

y

n n

x

n n

y

n

x

n

y

n

x

n

n x y x y x x y y

∆

− +

∆

− = ∆

∆

+

∆

∆ = ∂

∂

+

∂

∂ = ∇ • = 2 1 2 1 κ 2 . (45) 

The terms ∆x and ∆y are computed from the cell center locations as shown in Figure 3.40. The 

curvature at vertices (x3, y3), (x4, y4) is obtained in a similar manner. The face center curvature is 

calculated by taking the average of the curvature at nearby vertices. 

κ2

(x1,y1) 

(x2,y2)

(x3,y3)

x y 

n1 n2

(xc2,yc2) (xc3,yc3)

κ3

κ4

(xc4,yc4)

(x4,y4)

(x,y) : node position 

(xc,yc) : cell position 

∆x 

∆y 

Figure 3.40 Diagram of the algorithm for calculating surface curvature. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 59 

3.5.5 Parametric Study of Grid Sensitivity and Time Accuracy 

In Task 1 to Task 3, it was observed that the CPCFD simulation results are sensitive to grid 

resolution and time step size. In Task 4, a set of carefully chosen parametric cases have been 

performed to investigate the grid sensitivity and time accuracy of CPCFD before tackling the 

three required sub-tasks. These parametric cases are summarized in Table 3.12. 

Table 3.12 Parametric Cases for Studying Grid Sensitivity and Time Accuracy. 

Case Time Step Grid CPU time (hr) Surface tension effects 

T7 ∆tmax=0.1; CFL=0.5 400X30 8.5 Not included 

T18 ∆tmax=0.045; CFL=0.5 400X30 15.0 Not included 

T8 ∆tmax=0.02; CFL=0.5 400X30 15.0 Not included 

T13 ∆tmax=0.01; CFL=0.5 400X30 34.0 Not included 

T19 ∆tmax=0.02; CFL=0.5 533X30 15.0 Not included 

T14 ∆tmax=0.02; CFL=0.5 600X30 17.0 Not included 

T15 ∆tmax=0.02; CFL=0.5 800X30 21.0 Not included 

T9 ∆tmax=0.02; CFL=0.5 400X30 21.0 Included 

T16 ∆tmax=0.02; CFL=0.5 600X30 36.0 Included 

T17 ∆tmax=0.02; CFL=0.5 800X30 48.0 Included 

All cases were performed for the PP702N resin. The simulated time period is 250 seconds and all 

cases were run on a Dell Precision 450 Desktop PC with a Xeon 3 GHz CPU. Note that the timestep size is determined by two factors, namely the maximum time-step size ∆tmax and the CFL 

number. The grid refinement is only performed in the horizontal direction. 

The melt flow rate is a constant at 0.14 gs-1 , which is an estimate of the mean mass flow rate 

based on the simulation results from Task 3. In Case 3R (see Figure 3.34), the predicted mean 

total mass loss rate is 0.27gs-1, the fraction of melt flow occurring as melt flow is 0.74. 

Considering the object modeled in Task 3 has a width (W) of 10 cm and the object modeled in 

Task 4 has a with of 7 cm, the melt flow rate used in Task 4 can be calculated as 

0.27 × 0.74× 7 /10 = 0.14 gs-1. 

The mass sources are distributed in the 10 cells at the center of the top surface in the 400X30 

grid and the 20 center cells at the center of the top surface in the 800X30 grid. A 2.5 degree tilt 

angle is applied to the catch plate by adjusting the direction of the gravity force. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 60 

Figure 3.41 shows the impact of the time-step size on the free surface shape. The “ringing” 

effects observed at the two edges of the melt pool are non-physical; these are caused by 

numerical errors in the solution algorithm. It can be seen that, due to the tilted catch plate, the 

right edge of the melt pool spreads faster than the left edge. As the time-step gets smaller, the 

spreading velocity of the melt pool increases. 

To quantify the impact of the time-step size on the calculated melt pool behavior, the mean 

velocity of the right edge of the melt is plotted against the square of the maximum time-step size 

in Figure 3.42. The cases plotted are for cases T7, T18, T8 and T13. The mean velocity of the 

edge is obtained by dividing the moving distance of the edge by the time difference. It can be 

seen that the mean velocity at 

250 seconds increases linearly 

as the square of the maximum 

time-step size decreases, 

which indicates that the 

numerical scheme is second 

order accurate in time for this 

Task 4 problem. The timestep size used in Case T7 is 

too large and causes 

significant errors in predicting 

the melt pool spreading rate. 

Note that surface tension 

effects were not included in 

the four test cases. 

T7 T8

T8 T7

T = 50 sec. 

T = 125 sec. 

T8 T7 

T13 

T13 

T13 

T = 250 sec. 

Figure 3.41 Free surface shapes at selected times for different time-step sizes. 

2.80E-04

2.85E-04

2.90E-04

2.95E-04

3.00E-04

3.05E-04

3.10E-04

3.15E-04

3.20E-04

0 0.002 0.004 0.006 0.008 0.01 0.012

Squre of Maximum Time Step (∆tmax2

)

Mean Velocity of Right Front at 250 sec. (m/s)

Figure 3.42 Mean spreading velocity as a function of the 

square of the maximum time-step size. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 61 

Figure 3.43 and Figure 3.44 show the free surface shapes for simulations performed with 

different grid resolutions in the horizontal direction. The surface tension effects are included in 

the three cases shown in Figure 3.44 but not included in the three cases shown in Figure 3.43. 

Note that the time-step size criteria for all six cases are the same. The first observation from the 

two figures is that the surface tension effects make the shape at the edge of the melt pool 

smoother. The second observation is that the “ringing” phenomenon becomes more severe as the 

grid resolution increases. The “ringing” effects are similar to that of the well-known Gibbs 

phenomenon for discontinuous functions [14]. 

T8 

T15 

T = 100 sec. 

T8 

T15 

T = 200 sec. 

T14 

T14 

Figure 3.43 Free surface shapes at selected times using different grid. No surface tension 

effects included. 

T9 

T17

T = 50 sec. 

T9 

T17

T = 130 sec. 

T16

T16 

Figure 3.44 Free surface shapes at selected times using different grid. Surface 

tension effects included. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 62 

To quantify the impacts of the grid resolution on the predicted melt pool behavior, plotted in 

Figure 3.45 is the mean velocity of the right edge of the melt versus the square of the normalized 

grid size. The normalized grid size is defined as the horizontal grid cell size divided by the 

horizontal grid cell size for the 400 × 30 grid (i.e., 0.6mm). The data points for surface tension 

not included correspond to simulations T7, T19, T`14 and T15 and the data points for surface 

tension included correspond to simulations T9, T16 and T17. From the figure it can be seen that 

a linear relationship exists between the data points, indicating that the numerical scheme used in 

CPCFD is second order accurate in space. 

Figure 3.45 Mean spreading velocity as a function of the square of the normalized grid size. 

1.0E-04

1.5E-04

2.0E-04

2.5E-04

3.0E-04

3.5E-04

0 0.2 0.4 0.6 0.8 1 1.2

Mean Right Front Velocity at 130 sec. 

S urface Tension E ffects not Included

S urface Tension E ffects Included

“Accurate solution” by Extrapolation 

Square of Normalized Grid Size ( (∆X / ∆ X 400x30) ) 2
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 63 

3.5.6 Task 4A, 4B and 4C Results and Discussions 

In Tasks 4A-4C, simulations are performed for two polymer resins with different catch plate tilt 

angles. Three cases have been performed and are summarized in Table 3.13. A 400 × 30 grid is 

used in all cases, where the time dependent mass source is distributed among the ten center cells 

near the free surface. It was found that the small melt flow rate onto the catch plate at the 

beginning (i.e., < 100 sec.) caused the case to crash because the melt temperature became too 

low due to heat loss. To overcome this difficulty, the melt flow on the catch surface from 0 

seconds to 200 seconds is modeled using a constant melt flow rate in all cases to create an initial 

amount of mass on the catch plate. The simulations are performed in the following steps: 

1. Compute the melt flow using a large constant melt flow rate (i.e., 0.14gs-1) for 30 seconds. 

This will result in about 5 grams of melt on the catch plate. 

2. Stop simulation. 

3. Restart simulation using the time dependent melt flow rate profile, starting at 200 seconds. 

4. Continue simulation until melt front approaches end of plate. 

Table 3.13 Task 4 Case Summary. 

Case Resin Melt 

Flow 

Rate 

Time Step 

Size 

Tilt 

Angle 

Simulated Time 

Period 

Wall Clock Time 

4A PP702N From 

Case 3R 

∆tmax=0.02; 

CFL=0.5 

2.5° 200 s – 430 s ~19 hrs 

4B PP6523 From 

Case 3I 

∆tmax=0.02; 

CFL=0.5 

2.5° 700 s – 1350 s ~36 hrs 

4C PP702N From 

Case 3R 

∆tmax=0.01; 

CFL=0.5 

0.0° 200 s – 540 s ~40 hrs 

The time dependent melt flow rates are obtained from the Task 3 modeling results. For Case 4A 

and 4C, the melt flow rate is computed from the mass loss history predicted by Case 3R (see 

Figure 3.46); additional details on the Case 3R results can be found in Section 3.4.3 and Figure 

3.34. The melt flow rate used in Case 4B (see Figure 3.47) is calculated from the Case 3I results, 

which can be found in Section 3.4.4 and in Appendix A. 

In Figure 3.46 and Figure 3.47, the solid lines represent the total mass loss from the polymer 

slab. The dashed lines are the mass loss occurring as melt flow. The difference between the solid 

line and dashed line is the mass loss due to gasification. The melt flow rate is calculated as the 

first derivative of the dashed line and is stored in a look-up table, which is read-in as input data 

by CPCFD at the beginning of the simulation. During the simulation, the flow rate is tabulated 

from the table given at a certain point in time. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 64 

. 

0

0.05

0.1

0.15

0.2

0.25

0.3

0 100 200 300 400 500 600 700 800 900

Time (s)

Melt flow rate (g/s)

0

50

100

150

200

250

300

0 200 400 600 800

Time (s)

Mass Loss (g)

Case 3R 

Figure 3.46 Melt flow rate used in Case 4A and Case 4C.

0

50

100

150

200

250

300

0 500 1000 1500 2000

Time Exposed (s)

Mass Loss (g)

Case 3I 

0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

0 400 800 1200 1600

Time (s)

Melt Flow Rate (g/s)

Figure 3.47 Melt flow rate used in Case 4B. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 65 

Figure 3.48 shows the temperature contours of the melt pool and the catch plate calculated in the 

three cases. For all cases, the temperature at the center of the melt pool near the free surface is 

higher than the other part of the melt pool because the hot melt from the polymer slab enters the 

melt pool at this location. Away from the center, the free surface becomes cool due to heat loss 

to the environment. The melt pool is heated by the catch plate where a temperature gradient in 

the vertical direction can be seen. Due to the tilted catch plate in Case 4A and 4B, the right front 

of the melt pool moves faster than the left front. In Case 4C, the two fronts of the melt pool move 

at the same pace and the shape of the melt pool as well as the temperature distribution are 

symmetric with respect to the center of the domain in the horizontal direction. The surfaces in all 

three cases are very smooth. 

T (K)

Case 4B 

850 second 

1000 second

1350 second

T (K)

Case 4A 

250 second 

325 second

400 second

Case 4C 

250 second 

325 second 

400 second

T (K)

Figure 3.48 Temperature contours of the melt pool and catch plate. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 66 

Comparing Case 4A and 4B in Figure 3.48, it can be seen that: 1) the overall height of the melt 

pool in Case 4A is larger than that in Case 4B; and 2) the right front in Case 4A moves faster 

than the right front in Case 4B. This is caused by the different viscosity-temperature 

relationships of PP702N and PP6523, and the different melt flow rates used in the two cases. 

Illustrated in Figure 3.49 is the time history of the melt front velocities for Case 4A, Case 4B and 

Case 4C. The plots show that the right melt front moves faster than the left melt front for Case 

4A and Case 4B. However, for Case 4C the two melt fronts move at the same speed. In the 

future, when the data becomes available, it will be possible to compare these plots to NIST 

experimental data. 

The “zigzag” appearance of the curves plotted in Figure 3.49 is a result of the melt front tracking 

algorithm used in CPCFD. In the algorithm, the position of the melt front is “marked” by the 

coordinate of the grid point immediately behind the front. The “marked” position of the melt 

front remains unchanged until the melt front moves a distance ∆x (i.e., one grid cell in the 

horizontal direction). The mean velocity reported in Figure 3.49 is obtained using the following 

formula: 

 

0

0

t t

x x

V t t

front −

− = (46) 

where xt and xt0 are the positions of the front at the current time and at the beginning of the 

simulation, respectively. Because xt is not updated continuously (i.e., not updated at every time 

step), the velocity vs. time curves appear squiggly. 

Figure 3.49 Mean velocity of melt front as a function of time. 

0.0E+00

1.0E-04

2.0E-04

3.0E-04

4.0E-04

5.0E-04

200.00 250.00 300.00 350.00 400.00 450.00

Time (s)

Mean Velocity of Melt Front (m/s)

4A, Right Front

4A, Left Front

0.0E+00

2.0E-05

4.0E-05

6.0E-05

8.0E-05

1.0E-04

1.2E-04

1.4E-04

700 800 900 1000 1100 1200 1300 1400

Time (s)

Mean Velocity of Moving Fronts (m/s)

Right Front

Left Front

0.0E+00

5.0E-05

1.0E-04

1.5E-04

2.0E-04

2.5E-04

3.0E-04

3.5E-04

200 300 400 500 600

Time (s)

Mean Velocity of Moving Fronts (m/s)

Right Front

Left Front Case 4A Case 4B Case 4C
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 67 

4. CONCLUSIONS 

In this report is described the work effort by REI to develop, demonstrate and deliver to NIST a 

condensed phase CFD based tool to model the processes of melting, flow and gasification of 

thermoplastic materials exposed to a high heat flux. Potential applications of the tool include 

investigating the behavior of polymer materials commonly used in personal computers and 

computer monitors if exposed to an intense heat flux, such as occurs during a fire 

The model delivered to NIST is based on a time dependent (time varying) grid CFD method. 

• The model is written in FORTRAN 90 in an object-oriented form. A 3D, finite volume, 

multi-block body-fitted time dependent (time varying) grid formulation is used to solve 

the unsteady Navier Stokes equations. The time integration, spatial discretization and 

overall solution procedure are based on standard CFD methods from the literature. A 

multi-grid method is used to accelerate convergence at each time step. 

• Sub-models are included to describe the temperature dependent viscosity relationship and 

in-depth gasification and absorption of thermoplastic materials, free surface flows and 

surface tension. NIST data is used for key material properties of the thermoplastic 

materials of interest. 

• A variety of boundary conditions can be used for the velocity field (no-slip, free-slip) and 

heat transfer to the object (adiabatic, heat loss, specified heat flux). 

• Model outputs include the time dependent velocity, temperature and position 

(displacement) at points in the thermoplastic body which can be imported to standard 

CFD visualization packages. Additional outputs include the time history of the mass loss 

rate and heat fluxes. 

• The accuracy and capabilities of the modeling tool are demonstrated on a series of test 

cases of increasing complexity. The test cases include grid sensitivity studies, adding heat 

loss boundary conditions, simulations for two thermoplastic materials (PP702N, PP6523), 

different heat flux scenarios and test problem configurations. 

Comparisons of model results to NIST experimental data indicate discrepancies between the 

model and experimental results, particularly for the rate of mass loss. To match NIST data for 

mass loss rate large changes were required to the parameters originally used in the in-depth 

absorption model and kinetic rate parameters in the in-depth gasification model. In addition, for 

simulations in which the free surface of the melt flow is parallel to the direction of gravity a grid 

smoothing operation needed to be applied to the free surface to control grid skewness that would 

lead to simulation divergence. For simulations that did not include models for all key physics 

(e.g., no in-depth absorption, no in-depth gasification) the simulations would develop a large 

deformation ("belly") at the free surface. For carefully selected model parameters and 

procedures, comparisons of the model results and NIST data show favorable agreement. 

However, the research effort was un-able to provide satisfactory explanations for the need to 

significantly increase the model parameters for the in-depth absorption and in-depth gasification 

models or the presence of the "belly" at the free surface of the thermoplastic object. These issues 

remain un-resolved though some efforts continue at NIST to address these issues. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 68 

The current work effort has highlighted the importance of having accurate models to describe indepth absorption, in-depth gasification and the highly non-linear viscosity-temperature 

relationship that exists for thermoplastic objects. Hopefully, future experiments can provide the 

necessary data to develop more accurate models for these items. 

The source code, documentation, test problems (all needed files) and presentation material have 

been delivered to NIST. 

 

In the future, it may be possible to couple models such as the condensed phase CFD code to the 

NIST FDS code, a CFD model that solves for gas phase transport and combustion. The 

combined tools would enable high-fidelity simulations for fire spread scenarios in office 

environments populated with desktop/deskside personal computers. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 69 

5. LITERATURE REFERENCES 

1. Butler, K. M.; Ohlemiller, T. J.; Linteris, G. T., “Progress Report on Numerical Modeling of 

Experimental Polymer Melt Flow Behavior,” Interflam 2004. International Interflam 

Conference, 10th Proceedings, Volume 2, July 5-7, 2004. Edinburgh, Scotland, Interscience 

Communications Ltd., London, England, 937-948 pp, 2004. 

http://fire.nist.gov/bfrlpubs/fire04/PDF/f04032.pdf

2. Ohlemiller, T. J.; Shields, J. R.; Butler, K. M.; Collins, B.; Seck, M., “Exploring the Role of 

Polymer Melt Viscosity in Melt Flow and Flammability Behavior,” New Developments and 

Key Market Trends in Flame Retardancy. Fall Conference. Proceedings. Fire Retardant 

Chemicals Association. October 15-18, 2000, Ponte Vedra, FL, Fire Retardant Chemicals 

Assoc., Lancaster, PA, 1-28 pp, 2000. http://fire.nist.gov/bfrlpubs/fire00/PDF/f00128.pdf

3. Ohlemiller, T. J.; Butler, K. M., “Influence of Polymer Melt Behavior on Flammability,” 

U.S./Japan Government Cooperative Program on Natural Resources (UJNR). Fire Research 

and Safety. 15th Joint Panel Meeting. Volume 1. Proceedings. March 1-7, 2000, San 

Antonio, TX, Bryner, S. L., Editor(s), 81-88 pp, 2000. 

http://fire.nist.gov/bfrlpubs/fire00/PDF/f00129.pdf 

4. Yan, Z., “CFD Fire Simulation and Its Recent Development”, Dept. of Fire Safety 

Engineering, Lund University, Sweden. Chapter to be included in upcoming Developments in 

Heat Transfer Series published by WIT press.

5. Ferziger, J.H., and Peric, M., “Computational Methods For Fluid Dynamics”, Springer, 3rd

Edition, 2002. 

6. Patankar, S.V., “Numerical heat transfer and fluid flow”. McGraw-Hill, New York, 1980 

7. J.U. Brackbill, D.B. Kothe, & C. Zemach, “A continuum method for modeling surface 

tension”, J. Comput. Phy., 100, pp335-354, 1992. 

8. Tang, Q., “MEPDF Method for Modeling Turbulent Combustion Using Detailed Chemistry”, NSF 

SBIR Phase I Final Report, DMI-0441833, July, 2005. 

9. W. Zhao, C. Zhang, S. H. Frankel, & L. Mongeau. Computational Aeroacoustics of 

phonation, I: Numerical methods, acoustic analogy validation, and effects of glottal 

geometry. J. Acoust. Soc. Am., 112(5), pp2134-2146, 2002. 

10. W. Zhao. A numerical investigation of sound radiation from subsonic jets with application to 

human speech production, PhD thesis, Purdue University, West Lafayette, 2000. 

11. Pierre-Gilles de Gennes, Françoise Brochard-Wyart, David Quéré (2002). Capillary and 

Wetting Phenomena -- Drops, Bubbles, Pearls, Waves. Springer. 

12. W.A. Sirignano and I. Glassman, Combustion Sci. and Tech., Vol. 1, pp. 307-312, 1970. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects 70 

13. Thomas Ohlemiller, Kathy Butler, NIST, private communication. 

14. J.W., Gibbs, "Fourier Series". Nature 59, 200, 1898 and 606, 1899. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects A - 1 

APPENDIX A 

Detailed PowerPoint Summary of Modeling Results 
Task 0
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 0: port CPCFD to Fortran 90

(10/01/05 ~ 10/25/05)
0

0.010.020.030.040.050.06

0 500 1000 1500 2000

Time (s)

H (x=0) (m)

C++Fortran90

Object height as a function of time 
Temperatur

Object Temperature Contour 

time = 1000 s

C++ Code

Fortran90 Code
Velocity Vector Field 

|V| (m/s) time = 1000 s

C++ Code

Fortran90 Code
Summary

During this performance period, we haveconverted the C++ version CPCFD code to a Fortran 90 version compared the performance of the two codes on the triangle test 

 the two codes give exactly the same results (within the round-off 

precision floating point)

corrected several coding errors in the original C++ codeidentified several areas that need improvement or modification (s

work)

 these improvements and modifications are necessary for the melt-d

in task 1 and other future work
Future Work

Test numerical scheme for implicit updating free surface position

 Currently using explicit scheme that requires small time step to avo

Test generalized algorithms for tracking free surface position and 

grids

 Current algorithm is specifically designed for the triangle test proble

Test methods for assigning stress free boundary condition to free 

 Currently, free surface boundary is treated as outlet boundary (zero 

is not stress free (though zero mass flow rate across the free surfacguaranteed).

 The current method may cause problem in cases where surface ten

important or if the effects of surface tension are required to study in 

2D condensed phase melt-drip model with steady heat flux 
Task 1
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 1: 2D condensed phase melt-drip mo

with steady heat flux 

(1/13/06)
Status Summary

 Investigated “Model 2” (see Progress Report – Dec, 2005)

3⁄4 Solution diverged at about t = 230 s

 Investigated a simplified problem

3⁄4 Replaced bottom no-slip wall BC with an outlet boundary3⁄4 Used rectilinear grid 3⁄4 Simulation ran successfully

 Identified deficiencies in code and model that lead to insta

(grid skewness may not be the only reason!)

9 fluctuations in free surface shape9 time step too large9 localized over-heating of polymer in very thin regions

 Enhanced robustness of solver

9 a filtering function to suppress fluctuation of free sur9 improved time step size selection9 special treatment to avoid over-heating in very thin r
Recall “Model 2” and problems enco

Velocity vectors at the corner of the modelSolution diverged
A Simplified Problem

Left face of sampleConstant heat flux 

q

surface of flowink 

∂T/∂s = - q

0, wher

direction normal to thRight face of sampu = 0 (no penetratiov = 0 (no-slip)∂T/∂x = 0 (adiabatic

Top face of samplev = 0 (no penetratiou = 0 (no-slip)∂T/∂y = 0 (adiabatic

Bottom face of samZero gradient (outletAll other parametebefore
Rectilinear Grid and Solution 

U velocity

T = 710 s

V velocity Viscosity Temperature
Example – fluctuations of moving bo

T=140s T=260s

Use filterinsuppress t
Example - local over heating

Temperature

• Temperature can reach mo

than 2000 K in regions whepolymer is extremely thin

• Leads to instability• Results in very small time s• Special treatment:

when material thickness is 

than a small value, set the flux to zero 
Movies – polymer temperature

Temperature (K) Te

Without special treatmentTo avoid local over-heating

With special treatmeTo avoid local over-h
Mass Loss vs. Time

0

20406080

100

0 200 400 600 800 1000 1200

Time (s)

Mass Loss Percentage 

25X25020X20010X100
Velocity Profile at Outlet

-0.006-0.004-0.002

0.000

0 0.01 0.02

X (m)

V velocity (m/s) t = 300 s t = 500 s t = 700 s

t= 900 st = 1100 s

Non-zero velocity near free surface 

~zero flow 

Results for 10x100 grid
Next Steps

Comparison to NIST data

• NIST to provide data for mass loss vs. time • NIST experimental data• NIST VOF numerical model results (if available)

Use improved solver to re-visit original problem (with no-slip BC) u

single block rectilinear grid

If needed, switch to multi-block model

3⁄4 current work will be very useful to guide multi-block simulation
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 1: 2D condensed phase melt-drip mo

with steady heat flux 

(1/20/06)
Status Summary

 Investigated 2D problem with no-slip wall BC at bottom of slab (ho

3⁄4 used improved code (see 11/13/06 report)3⁄4 modified re-meshing strategy near bottom wall boundary3⁄4 modified implementation of the outlet boundary condition3⁄4 100x10 grid used to investigate feasibility of model3⁄4 model with biased grid (i.e., fine grid near slab surface) in pro
Problem

Left face of sampleConstant heat flux 

q

to surface of flowing k 

∂T/∂s = - q

0, whe

the direction normsurfaceRight face of sampu = 0 (no penetratiov = 0 (no-slip)∂T/∂x = 0 (adiabatic

Top and Bottom samplev = 0 (no penetratiou = 0 (no-slip)∂T/∂y = 0 (adiabatic
Rectilinear Grid and Solution 

U velocity

T = 500 s

V velocity Viscosity Tempe

No-slip Wall Outlet Boundary

Initial Grid
Movies – polymer temperature

Temperature Temp

With no-slip bottom wall Without no-slip bottom w
Mass Loss vs. Time

0

102030405060708090

100

0 200 400 600 800 1000 1200

Time (s)

Mass Loss Percentage
Ave. Surface Temperature vs. Time

0

50

100150200250300350400450500

0 500 1000 1500 2000

Time (s)

Average Surface Temperature (°C)
Comparisons

Slip and No Slip Models

3⁄4 Flow begins at about t ~ 200 sec3⁄4 ~95% Mass Loss at t ~ 1200 sec

 No Slip – time to achieve 95% Mass Loss less for finer gr

 Ave. Surface Temperature : 250 – 450 °C during mass loss t

NIST Test Data

Flow starts at t ~ 200 sec~95% Mass Loss at t ~ 2300 sec Ave. Mass loss rate ~ 0.12 g/sec Surface Temperature ~ 350 °C

Physics Not Included in current model

9 radiative, convective heat transfer at surface9 gasification9 vaporization of polymer 
Next Steps

Complete grid sensitivity study for current polymerComparison to NIST data

• Compare mass loss or surface temperature vs. time ? Other ?

3⁄4 NIST experimental data = what is thickness of test piece 

Repeat tests using improved solver for second polymer

3⁄4 Need NIST data for material properties

Conf. Call by Jan. 30Complete Task 1 Jan. 31
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 1: 2D condensed phase melt-drip mo

with steady heat flux 

(1/26/06)
Modeled Problem

Left face of sampleConstant heat flux 

q

to surface of flowing k 

∂T/∂s = - q

0, whe

the direction normsurface (

q0 = 20000 

Right face of sampu = 0 (no penetratiov = 0 (no-slip)∂T/∂x = 0 (adiabatic

Top and Bottom samplev = 0 (no penetratiou = 0 (no-slip)∂T/∂y = 0 (adiabatic
Computation Domain and Grids

No-s

• 1• B

SD

• U

SD

No-slip wall

• 10X100 cells• Uniformed Grid 

Spacing in Horizon Direction

• Uniformed Grid 

Spacing in Vertical Direction

• 16X100 cells• Biased Grid 

Spacing in Horizon Direction

• Uniformed Grid 

Spacing in Vertical Direction

No-slip wall

* Grid lines in the plots connect cell centers
Residual History

-6-5-4-3-2-1

012345

0 50 100 150 200 250

Number of Iteration

Residual

MassU VelocityV VelocityTemperature

A typical residual history in one time step

⎟⎟⎟⎠⎞ ⎜⎜⎜⎝⎛ ≡ ∑∑ Calculatoin domain Y Calculation domain 0 10 log ψ ψ Y Y R R

Residual definition

Where RY ≡AEYE +AWYW +ANYN +ASYS +ATYT +ABYB +SU −AP

0

And ψ Y = APYP

Y is the solved quantity

• D

itewsobesirefrite
Computational Cost

• All Cases were run on a Dell Precision 450 Desktop PC with a Xeon 3GHz CPU

14X200 12 hours ~2 MB14X100 6 hours ~1 MB10X100 4 hours < 1MB

Case Wall Clock Time Memory Requirement
T=10s T=50s T=100s T=150s T=200s T=250s T=300s T=350s T=400s TFree Surface Location and Temperature Field Against Time

0% mass loss 9
V velocity (

T=10s T=200s T=250s T=300s T=350s T=400s T=450s

Free Surface Location and Velocity Field Against Time

• Negative vel

downward

• Melted mate

layer and flo

• “Gold” colo

small velocifree surface 

Grid: 14X
Velocity Vectors at Time = 300 s

Grid: 14X200
Temperature (K)

Grid: 14X200

Movie: Free Surface Location and Temperature Contour

• Free 

(wavoutleupwa
0

102030405060708090

100

0 100 200 300 400 500

Time (s)

Percentage of Mass Loss

14X20014X10010X100

Mass Loss as Function of Time

14X200 21.4* Computed from 200 s to 450 s

Experiment 0.12

14X100 25.4* Computed from 200 s to 410 s10X100 28.1* Computed from 250 s to 440 s, 

Grid Mass Loss Rate Comments

* per unit length (g/s m)
0.00E+00

2.00E-034.00E-036.00E-038.00E-031.00E-021.20E-021.40E-021.60E-021.80E-022.00E-02

0 100 200 300 400 500

Time (s)

Velocity Magnitute on Free Surface (m/s)

14X20014X10010X100

Area Weighted Average Velocity on Free Surface as Function of 
Area Weighted Average Temperature on Free Surface as Functi

14X200 347 °C

Experiment 340 °C

14X100 303 °C Averaged from 0 to 95% mass loss 10X100 295 °C

Surface Comments

Temperature

Grid

0

50

100150200250300350400450500

0 100 200 300 400 500

Time (s)

Temperature on Free Surface (°C)

14X20014X10010X100
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 1: 2D condensed phase melt-drip mo

with steady heat flux (PP702) 

(3/06/06)
Modeled Problem

Left face of sampleConstant heat flux 

q

to surface of flowing k 

∂T/∂s = - q

0, whe

the direction normsurface (

q0 = 20000 

Right face of sampu = 0 (no penetratiov = 0 (no-slip)∂T/∂x = 0 (adiabatic

Top and Bottom samplev = 0 (no penetratiou = 0 (no-slip)∂T/∂y = 0 (adiabatic
Solution Algorithm

1. Solve the momentum equations using the geometry defined by the current shape of 

surface and the prescribed pressure at it

2. Enforce local mass conservation in each CV by solving the pressure-correction equa

prescribed pressure boundary condition on the current free surface. Mass is conservglobally and in each CV, but the non-zero mass fluxes through the free surface may 

3. Correct the position of the free surface so that the volume defined by its corrected an

position (i.e., at previous time step) compensates the mass fluxes through the free s

obtained in the preceding step.

4. Re-mash the grids based on the corrected free surface position and update the trans

Jacobian matrix.

5. Return to step 1 and repeat until all equations and boundary conditions are satisfied

6. Advance to the next time step
Computation Domain and Grids

No-s

• 1• B

SD

• U

SD

• G

0h

No-slip wall

• 10X100 cells• Uniformed Grid 

Spacing in Horizon Direction

• Uniformed Grid 

Spacing in Vertical Direction

• Grid size less than 

2.5 mm in horizontal direction

• 14X200 cells• Biased Grid 

Spacing in Horizon Direction

• Uniformed Grid 

Spacing in Vertical Direction

• Grid size less than 

1.0 mm in horizontal direction

No-slip wall

* Grid lines in the plots connect cell centers
Residual History

-6-5-4-3-2-1

012345

0 50 100 150 200 250

Number of Iteration

Residual

MassU VelocityV VelocityTemperature

A typical residual history in one time step

⎟⎟⎟⎠⎞ ⎜⎜⎜⎝⎛ ≡ + ∑∑ Calculatoin domain Y Calculation domain 0 10 7 log ψ ψ Y Y R R

Residual definition

Where RY ≡AEYE +AWYW +ANYN +ASYS +ATYT +ABYB +SU −AP

0

And ψ Y = APYP

Y is the solved quantity

• D

itewsobesirefrite
Computational Cost

• All Cases were run on a Dell Precision 450 Desktop PC with a Xeon 3GHz CPU

19X200 ~100 hours ~2 MB14X200 ~ 30 hours ~1 MB

10X100 ~ 4 hours < 1MB

Case Wall Clock Time Memory Requirement
T=10s T=50s T=100s T=150s T=175s T=200s T=215s T=230s T=240s TFree Surface Location and Temperature Field vs. Time

0% mass loss 8
Velocity Ma

T=10s T=200s T=210s T=220s T=230s T=240s T=250s

Free Surface Location and Velocity Magnitude vs. Time

Grid: 19X

• Melted mate

layer and flo
Velocity Vectors at Time = 300 s

Grid: 19X200
Temperature (K)

Grid: 19X200

Movie: Free Surface Location and Temperature Contour
Volume/Mass Change as Function of Time

0

0.00010.00020.00030.00040.00050.00060.0007

0 50 100 150 200 250 300

time (s)

volume (m3)

volume_remainvolume_outvolume_totalinitial volume 

0

102030405060708090

100

0 50 100 150 200 250 300

time (s)

percentage

mass losscumulative error in total mass

• Cumu

mass definenormabetwewithin the suof the mass 

• Cumu

than 1reachiand is 95% mincreaerror in conseof the to the skewn
Mass Loss as Function of Time

* per unit length (g/s m)

0

102030405060708090

100

0 100 200 300 400 500

time (s)

Mass loss percentage

19X20014X20010X100

• The mass loss ra

m) in the case wicells, while the nexperiment is 0.1

• The three cases 

grid convergencethe grids beyond computationally pthis point to invesissue.
0

0.05

0.1

0.15

0.2

0.25

0 0.02 0.04 0.06 0.08 0.1 0.1

Velocity Magnitude (m/s)

Elevation (m)

200 sec (~10% mass loss)230 sec (~40% mass loss)240 sec (~60% mass loss)

Velocity Magnitude Profiles on Free Surface at Different Time
0

0.05

0.1

0.15

0.2

0.25

600 650 700 75

Temperature (K)

Elevation (m)

200 sec (~10% mass loss)230 sec (~40% mass loss)240 sec (~60% mass loss)

Temperature Profiles on Free Surface at Different Times

• Average free surface temperature observed in experiment was around 6

• The computation model does not include processes such as radiant and 

convective heat loss, gasification, etc., which may be the cause of the htemperature calculated on the free surface. 
A Computational Model For F

Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 1: 2D condensed phase melt-drip mo

with steady heat flux

(4/03/06)
Outline

• Problem Definition• Solution Algorithm• Three Computational Grids• Example Convergence History• Computational Cost• Example Results

– Free surface location, temp., velocity - selected times– Movie of free surface location and temperature– Mass loss, surface temperature, surface velocity

• Summary
Modeled Problem

Left face of sampleConstant heat flux 

q

to surface of flowing k 

∂T/∂s = - q

0, whe

the direction normsurface (

q0 = 20000 

Right face of sampu = 0 (no penetratiov = 0 (no-slip)∂T/∂x = 0 (adiabatic

Top and Bottom samplev = 0 (no penetratiou = 0 (no-slip)∂T/∂y = 0 (adiabatic
Initial conditions:T0 = 25 o C

u = 0v = 0p = 105 Pa

Problem parameters:q0 = 20,000 W/m2g = 9.8 m/s2hs = 0.25 m

ws = 0.025 m

Material properties: All parameters are constant except for viscosity.ρ = 900 kg/m3cp = 2400 J/kg-K

k = 0.25 W/m-K

Viscosity:

( )( ) ( )

( )

( )

( ) 

f

f T

f T

T f

⎪⎪⎪⎩⎪⎪⎪⎨⎧ − − + =

425

10 200 / 200 25 200

10

2

2

1

1

6

6

μ

( ) ( 

4 

f1 T 10^ 14.48 0.13858T 5.5960 10 

− 

= − + × 

( ) ( 

f 2 T 10^ 53.19 0.2542T 2.9879 

= − + 

Parameter Values – PP702 Resin
Initial conditions:T0 = 25 o C

u = 0v = 0p = 105 Pa

Problem parameters:q0 = 20,000 W/m2g = 9.8 m/s2hs = 0.25 m

ws = 0.025 m

Material properties: All parameters are constant except for viscosity.ρ = 900 kg/m3cp = 2400 J/kg-K

k = 0.25 W/m-K

Viscosity:

Parameter Values – PP6523 Resin

( )( ) ( )

( )

( )

( ) 

f

f T

f T

T f

⎪⎪⎪⎩⎪⎪⎪⎨⎧ − − + =

3

2

430

10 220 / 220 25 220

10

2

2

1

1

6

6

μ

( ) ( 

3 2 

f1 T 10^ 85.843 1.3870T 7.979 10 T 2.0286 

− 

= − + − × + × 

( ) ( 

f 2 T 10^ 100.90 0.4771T 5.533 

= − + 
Solution Algorithm

1. Solve the momentum equations using the geometry defined by the current shape of 

surface and the prescribed pressure at it

2. Enforce local mass conservation in each CV by solving the pressure-correction equa

prescribed pressure boundary condition on the current free surface. Mass is conservglobally and in each CV, but non-zero mass fluxes through the free surface may res

3. Correct the position of the free surface so that the volume defined by its corrected an

position (i.e., at previous time step) compensates the mass fluxes through the free s

obtained in the preceding step.

4. Re-mash the grids based on the corrected free surface position and update the trans

Jacobian matrix. 

• Grid smoothing to avoid gridline entanglement

5. Return to step 1 and repeat until all equations and boundary conditions are satisfied6. Advance to the next time step

Discretization Schemes

 Second-order central difference scheme used for spatial discretization of momentum  Second-order Total Variation Diminishing (TVD) scheme used for spatial discretizati

equation

 Second-order (three level) scheme used for time integration
Heat Sink Term for Energy Equation

A heat sink term in Arrhenius form is added to the energy equation to allev

overheating

• Overheating results in non-physical temperatures along the free surface due to 

physical models, which in-turn creates sharp changes in viscosity (e.g., one ord

change in viscosity due to several degrees difference in temperature) along the 

The material with low viscosity tends to roll-over the material with high viscosity

be handled by the current algorithm.

• The material roll-over may not cause problem for the VOF method.

• The overheating problem becomes significant when grid size becomes smaller n

surface due to the large temperature gradient near the surface

0.00E+005.00E+071.00E+081.50E+082.00E+082.50E+083.00E+083.50E+084.00E+084.50E+08

0 100 200 300 400 500 600 700 800

T (K)

q

• The heat sink term in the energy equation 

to mimic the gasification heat loss term. Hmass loss is involved in the model currentsink term can be expressed as

• The left figure shows the model used in th

presented in this report. In the model, the negligible when temperature is below 700 becomes significant when temperature exWe found this model works fine for PP702 which have similar melting temperatures aDifferent model must be used for PP23 bemuch lower melting temperature at 460 K

q

=

HρBexp

• Melting temperature can be found from the viscosity-temperature relationship wh

viscosity becomes constant when temperature exceeds the melting temperature
Computation Domain and Grids

* Grid lines in the plots connect cell centers

No-slip wall• 19X200 cells• Biased Grid Spacing in 

Horizontal Direction

9 smallest grid cell 

less than 0.5 mm

9 four grid cells at 

outlet boundary

9 Initial grid equivalent 

to 50X200 uniform grid in VOF

• Uniform Grid Spacing in 

Vertical Direction

5.e-4m

• 21• Bi

Ho

• Un

Ve

• 15X200 cells• Biased Grid Spacing in 

Horizontal Direction

9 smallest grid cell 

less than 0.625 mm

9 three grid cells at 

outlet boundary

9 Initial grid equivalent 

to 40X200 uniform grid in VOF

• Uniform Grid Spacing in 

Vertical Direction
Case Summary

PP23 10X100 (uniform) 1 2.5 10X100 FailedPP23 14X200 (biased) 2 1.25 20X200 Done PP23 21X200 (biased) 6 0.3125 80X200 5% mass loss PP23 19X200 (biased) 4 0.5 50X200 10% mass losPP702 21X200 (biased) 6 0.3125 80X200 DonePP702 19X200 (biased) 4 0.5 50X200 DonePP702 15X200 (biased) 3 0.625 40X200 DonePP702 14X200 (biased) 2 1.25 20X200 DonePP702 10X100 (uniform) 1 2.5 10X100 Done PP6523 21X200 (biased) 6 0.3125 80X200 DonePP6523 19X200 (biased) 4 0.5 50X200 DonePP6523 15X200 (biased) 3 0.625 40X200 DonePP6523 14X200 (biased) 2 1.25 20X200 Done

Status as of 03/29/06

Equiv. Uniform Grid in VOF

Initial Finest Grid Size (mm)

# of cells at outlet boundary

Polymer Grid
-6-5-4-3-2-1

012345

0 50 100 150 200 250

Number of Iteration

Residual

MassU VelocityV VelocityTemperature

Residual History

A typical residual history in one time step

⎟⎟⎟⎠⎞ ⎜⎜⎜⎝⎛ ≡ + ∑∑ Calculatoin domain Y Calculation domain 0 10 7 log ψ ψ Y Y R R

Residual definition

Where RY ≡AEYE +AWYW +ANYN +ASYS +ATYT +ABYB +SU −AP

0

And ψ Y = APYP

Y is the solved quantity

• P• D

itewsobesirefrite
Computational Cost

• All Cases were run on a Dell Precision 450 Desktop PC with a Xeon 

19X200 ~100 hours ~2 MB15X200 ~ 36 hours ~1 MB21X200 ~120 hours ~2 MB14X200 ~ 24 hours ~1 MB

10X100 ~ 4 hours < 1MB

Case Wall Clock Time Memory Requirement

PP702

• Run times are comparable between cases of PP702 and PP6523 wit
T=10s T=50s T=100s T=150s T=175s T=200s T=215s T=230s T=240s 

Free Surface Location and Temperature Field vs. Time – PP702
T=10s T=50s T=100s T=150s T=175s T=200s T=215s T=230s T=240s 

Free Surface Location and Temperature Field vs. Time – PP6523
Velocity Ma

T=10s T=200s T=210s T=220s T=230s T=240s T=250s

Free Surface Location and Velocity Magnitude vs. Time – 702

• Melted mate

layer and flo

Grid: 
Velocity Ma

T=10s T=200s T=210s T=220s T=230s T=240s T=250s

Free Surface Location and Velocity Magnitude vs. Time – PP652

• Melted mate

layer and flo

Grid: 
Example: Velocity Vectors at Time = 220 s (Zoom to 200% for be

PP7Grid

Velocity ~ except nea
Example: Velocity Vectors at Time = 235 s (Zoom to 200% for be

PP6Grid

• Velocity ~ except nea• Comparinflowing laythinner duviscosity atemperatur
Temperature (K)

Movie: Free Surface Location and Temperature Contour

PP702Grid: 19X200

Temperature 

PP6523Grid: 19X200
Volume/Mass Change as Function of Time

0

0.00010.00020.00030.00040.00050.00060.0007

0 50 100 150 200 250 300

time (s)

volume (m3)

volume_remainvolume_outvolume_totalinitial volume 

0

102030405060708090

100

0 50 100 150 200 250 300

time (s)

percentage

mass losscumulative error in total mass

• Cumu

mass definenormabetwewithin t=0) anmass and thdomai

• Cumu

than 1reachiand is 95% mincreaerror in conseof the to the skewnbottom 

PP7Gri
Volume/Mass Change as Function of Time

• Cumulative 

than 1% bereaching 82and is less 92% mass increasing error in maconservatioof the calcuto the signiskewness abottom lip 

PP6523Grid: 19

0

102030405060708090

100

0 100 200 300

Time (sec)

Percentage (%)

Mass LossCumulative Error in Mass
Mass Loss as Function of Time – PP702

21X200 62.4* Computed from 180 s to 260 s19X200 63.3* Computed from 180 s to 260 s

Experiment 0.12** Ohlemiller FAX 01-18-06

15X200 31.5* Computed from 210 s to 320 s14X200 28.1* Computed from 220 s to 330 s, 

“steady state” Comments

Mass Loss Rate

Grid

* per unit length (g/s m)** g/s

• The results are vegrid resolution near The resolution of thare not adequate to temperature gradie• The two cases witsimilar results in terhistory and “steady rate, which indicateto some extent.• The code chooses based on a constanthree cases. The timconstant in each cavery different in diffgeneral, the 21X20smallest time step serror also plays a roloss history differen• The predicted manot be directly comexperimental obseris a 2-D simulation; physical models are numerical model suradiation absorption

PP7

0

102030405060708090

100

0 50 100 150 200 250 300 350 400

Time (s)

Mass Loss Percentage (%)

PP702, 21X200PP702, 19X200PP702, 15X200PP702, 14X200
Mass Loss as Function of Time – PP6523

21X200 70.5* Computed from 200 s to 270 s19X200 70.7* Computed from 200 s to 270 s

Experiment N/A

15X200 59.1* Computed from 200 s to 290 s14X200 33.8* Computed from 250 s to 400 s, 

“steady state” Comments

Mass Loss Rate

Grid

* per unit length (g/s m)** g/s

• The case with coavery different meltincomparing to the ottook much longer tito flow and the masmuch smaller. • The differences acases with finer gridthan the differences the case with the coof the starting time ‘steady state’ mass loss history, which convergence to som

PP6

0

102030405060708090

100

0 50 100 150 200 250 300 350 400 450

Time (sec)

Percentage of Mass Loss (%)

PP6523, 21X200PP6523, 19X200PP6523, 15X200PP6523, 14X200
0

0.05

0.1

0.15

0.2

0.25

0 0.02 0.04 0.06 0.08 0.1 0.1

Velocity Magnitude (m/s)

Elevation (m)

200 sec (~10% mass loss)230 sec (~40% mass loss)240 sec (~60% mass loss)

Velocity Magnitude Profiles on Free Surface at Selected Times 
Velocity Magnitude Profiles on Free Surface at Selected Times 

0

0 .0 5

0 .1

0 .1 5

0 .2

0 .2 5

0 0 .0 5 0 .1 0 .1 5 

V e lo c ity M a g n itu d e (m /s )

Elevation (m)

2 1 5 s e c (~ 1 0 % m a s s lo s s )2 4 0 s e c (~ 4 0 % m a s s lo s s )2 6 0 s e c (~ 6 0 % m a s s lo s s )
0

0.05

0.1

0.15

0.2

0.25

600 650 700 75

Temperature (K)

Elevation (m)

200 sec (~10% mass loss)230 sec (~40% mass loss)240 sec (~60% mass loss)

Temperature Profiles on Free Surface at Selected Times

• Average free surface temperature observed in experiment (Ohlmeir FAX 

was ~613 K

• The computation model does not include processes such as radiant and 

heat loss, gasification, etc., which may be the cause of the higher tempe

calculated on the free surface. 
Temperature Profiles on Free Surface at Selected Times 

0

0 .0 5

0 .1

0 .1 5

0 .2

0 .2 5

550 600 650 700

T e m p e ra tu re (K )

Elevation (m)

2 1 5 s e c (~ 1 0 % m a s s lo s s )2 4 0 s e c (~ 4 0 % m a s s lo s s )2 6 0 s e c (~ 6 0 % m a s s lo s s )
Summary

• Completed Two Resins on Multiple Grids• Grid Sensitivity / Convergence

– 15x200 possibly too coarse– 19x200 and 21 x200 => about same mass loss rate

• Computational Cost highlighted 

– About same cost for both resins

• Model features highlighted

– Solution algorithm– Typical convergence history– Solutions at selected times– Velocity field resolution across melt boundary layer

• Example Results

– Free surface location, temp., velocity - selected times– Movie of free surface location and temperature– Mass loss, surface temperature, surface velocity
Task 2
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 2: 2D model with steady heat flux, ra

and convective heat losses, and in-dep

absorption of radiative heat flux

(5/31/06)
Radiative and convective heat loss, In-depth absorption models

Only impacts energy equationBoth radiative and convective heat losses are only applied to the free surfacA heat source term is added to the energy equation to account for in-depth aAbsorption coefficient is a function of the distance from the free surface

Absorption coefficient models

In-depth absorption coefficient data for one polymer resin provided by NIST

Email sent by Kathy to Bockelie on 2/22/2006Data applicable to PP23K in the HIFT apparatus

Three models constructed based on the NIST data

Cases studied

Computational gridCase summaryExample results: 

Free surface location, temp., velocityMass loss, surface temperature, surface velocity, and comparisons with data

Summary

Outline
Heat loss and In-depth absorption models

I0

qrad

I

Figure 1, Sketch of in-depth 

absorption and heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy equation source terms due to in-depth absorption, radiative and convecti

In-depth absorption: , assuming absorption process in x-direction oRadiative heat loss: , assuming radiative heat loss only on fConvective heat loss:Parameters: A – surface area; - volume; - emissivity (=1.0); - Stefan-Boltzmann Constan

T0 – ambient temperature (=298 K); hconv – convective heat transfer coefficient (= 

A heat sink term is added to the energy equation to prevent overheating (i.e., T > 750K) The models are implemented in subroutine fblck_bc_sca

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

σ
0

0.20.40.60.8

1

012345

x, Distance from Free Surface (cm)

Normalized Incident Heat Flux I/I

0Data provided by NIST

Model 1Model 2

Absorption Coefficient Models

Original data (symbols in Figure 2) is roughly applicable to the PP23K resin at the 

• Three incident heat flux distribution functions are used in the CFD model 

• Model 1 is differentiable, and all models are well defined over the entire thickness of the sample 

• Model 1: • Model 2:• Model 3: use the original fitted data

Figure 2, Incident heat flux distribution in the horizontal direction 

⎪⎩⎪⎨⎧

= ≤ <= ≤ <

− − +− +

I/I 10 , 0.3 x 1.3 (cm)I/I 10 , 0 x 0.3 (cm)

( 0.5154 1.473x 0.4339x )

0

( 5.17x 6.964x )

0

2

2

I/I exp( 3.843x 38.36x ), 0 x 3.0 (cm) 2 0 = − − ≤ ≤

⎪⎩⎪⎨⎧

= − − − ≤ <= − + ≤ <

I/I exp( 0.6814 0.7802x 0.009333x ), 0.3 x 5.0 (cm)

I/I exp( 5.17x 6.964x ), 0 x 0.3 (cm)

2

0

2

0
Computational Grid

Polymer sample dimension: 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Figure 3, Computational grid at t = 0 sec. 
Case Summary

Case 2D PP702 Included Included Model 2 27%, 480 sec ~ 2 h 

N/A~ 16 h~ 30 h~ 48 hEstimated CPU Time

N/A55%, 460 sec58%, 650 sec65%, 400 secMass Loss and Elapsed Time

Case 2C PP6523 Included Included Model 1 Case 2ECase 2BCase 2ACase

PP702 Included Included Model 3 PP702 Included Included Model 1 PP702 Included Not included N/A 

In-depth Absorption Coefficient Model

In-depth Absorption

Radiative and Convective Heat Losses

PolymerResin

All cases are run with a external Incident heat flux of 30KW/m2

* Manually stopped because calculation become extremely slow when mass loss exceeds 50%
Impact of In-depth Absorption on Melting Proces

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-

Time interval between movie frames: 5 secondsThe inclusion of in-depth absorption leads to significant difference in shap

Melting process slows down Upper free surface moves much slower towards the back wallMelting layer becomes thicker

Temperature contour
Impact of Different Resin on Melting Process 

Case 2B (PP702, heat loss + in-depth absorption) Case 2C (PP6523, heat loss + in-d

Time interval between movie frames: 5 secondsPP6523 has a much lower viscosity at high temperature, which leads to faprocess

Temperature contour
Impact of Different Absorption Coefficient Model

Case 2B (PP702, in-depth absorption model 1)

Temperature contour

Case 2D (PP702, in-depth absorpti

Time interval between movie frames: 5 secondsModel 2 (Case 2D) has a much larger optical thickness than that in Model 

deeper penetration of incident heat flux below the surfacethicker melting layer and slower melting processlarge amount of material bulge over the bottom lip and causes the calculation to 
Example: Velocity Vectors at t = 300 sec (zoom to 200% 

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-d

The addition of in-depth absorption leads to

slower melting process thicker melting layer below free surface

~50% mass loss ~15% mass loss
Mass Loss as Function of Time 

0

10203040506070

0 100 200 300 400 500 600 700 800

Time (s)

Mass Loss Percentage (%) Case 2A Case 2B Case 2C

Case 2D

Case 2D 14.8* Computed from 340 s to 480 s

Case 2C 18.0* Computed from 220 s to 440 s

0.15** Ohlemiller FAX 01-18-06 + Experiment

Case 2B 11.2* Computed from 220 s to 660 s

Case 2A 36.0* Computed from 200 s to 300 s, 

“steady state” Comments

Mass Loss Rate

Grid

* per unit length (g/s m) ** g/s + measured on a 2.5cmX25.0cm PP70
Example: Surface Temperature and Velocity Mag

0

0.020.040.060.08

0.1

400 450 500 550 600 650 700 750

Surface Temperature (K)

Elevation (m)

Case 2A, 280 secCase 2B, 440 secCase 2C, 350 sec

0

0.020.040.060.08

0.1

0 0.01 0.02 0.0

Surface Velocity Magnitud

Elevation (m)

Case 2ACase 2BCase 2C

35% mass l

35% mass loss

Surface temperature and velocity magnitude profiles at 35% mass 

Velocity magnitudes approach zero at upper no slip wallTemperature monotonically distributed along the surface in Case 2B and 2CInclusion of in-depth absorption decreases surface temperature and velocityOver-heating (T > 698K) occurs in Case 2A at upper portion of the free surface Case 2C has higher surface temperature and velocity than Case 2B due to low 

Measured surface temperature on a 2.5cmX25.0cm PP702 sample was 648 K (Ohlemiller FAX 01-18-06)
Summary

Studied Two Resins with Heat Losses and In-depth Absorption of Heat FluMost Cases were run to at least 50% mass loss

Calculation becomes extremely slow when mass loss exceeds 50%One case (Case 2D) diverged at 27% mass loss due to large optical thickness cOne case (Case 2E) is in progress using fitted absorption coefficient data

Computational Cost highlighted 

Less expensive than Task 1

Smaller computational domainLarger time step due to lower velocity caused by heat losses and in-depth absorption

Model features highlighted

Heat loss models and in-depth absorption modelsImpacts of different models on melting processesSolutions at selected timesVelocity field resolution across melt boundary layer

Example Results

Free surface location, temp., velocity - selected timesMovie of free surface location and temperatureMass loss, surface temperature, surface velocity, and comparison with experime
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 2: 2D model with steady heat flux, ra

and convective heat losses, and in-dep

absorption of radiative heat flux

(6/01/06)
Radiative and convective heat loss, In-depth absorption models

Only impacts energy equationBoth radiative and convective heat losses are only applied to the free surfacA heat source term is added to the energy equation to account for in-depth aAbsorption coefficient is a function of the distance from the free surface

Absorption coefficient models

In-depth absorption coefficient data for one polymer resin provided by NIST

Email sent by Kathy to Bockelie on 2/22/2006Data applicable to PP23K in the HIFT apparatus

Three models constructed based on the NIST data

Cases studied

Computational gridCase summaryExample results: 

Free surface location, temp., velocityMass loss, surface temperature, surface velocity, and comparisons with data

Summary

Outline
Heat loss and In-depth absorption models

I0

qrad

I

Figure 1, Sketch of in-depth 

absorption and heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy equation source terms due to in-depth absorption, radiative and convecti

In-depth absorption: , assuming absorption process in x-direction oRadiative heat loss: , assuming radiative heat loss only on fConvective heat loss:Parameters: A – surface area; - volume; - emissivity (=1.0); - Stefan-Boltzmann Constan

T0 – ambient temperature (=298 K); hconv – convective heat transfer coefficient (= 

A heat sink term is added to the energy equation to prevent overheating (i.e., T > 750K) The models are implemented in subroutine fblck_bc_sca

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

σ
Absorption Coefficient Models

Original data (symbols in Figure 2) is roughly applicable to the PP23K resin at the 

• Three incident heat flux distribution functions are used in the CFD model 

• Model 1 is differentiable, and all models are well defined over the entire thickness of the sample 

• Model 1: • Model 2:• Model 3: use the original fitted data

Figure 2, Incident heat flux distribution in the horizontal direction 

⎪⎩⎪⎨⎧

=

≤

<

=

≤

<

− − +− +

I/I 10 , 0.3

x 1.3 (cm)

I/I 10 , 0

x 0.3 (cm)

( 0.5154 1.473x 0.4339x )

0

( 5.17x 6.964x )

0

2

2

I/I exp( 3.843

x 38.36

x ),

0

x 3.0 (cm) 2

0 =

−

−

≤

≤

⎪⎩⎪⎨⎧

=

−

−

−

≤

<

=

−

+

≤

<

I/I exp(

0.6814

0.7802

x

0.009333

x ), 0.3

x 5.0 (cm)

I/I exp( 5.17

x

6.964

x ),

0

x 0.3 (cm)

2

0

2

0

0

0.20.40.60.8

1

012345

x, Distance from Free Surface (cm)

Normalized Incident Heat Flux I/I 0Data provided by NIST Model 1 Model 2 Model 3
Computational Grid

Polymer sample dimension: 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Figure 3, Computational grid at t = 0 sec. 
Case Summary

Case 2E PP702 Included Included Model 3 55%, 700 sec ~ 21 h 

Case 2D PP702 Included Included Model 2 27%, 480 sec ~ 2 h 

~ 7 h~ 16 h~ 30 h~ 48 hEstimated CPU Time

57%, 570 sec55%, 460 sec58%, 650 sec65%, 400 secMass Loss and Elapsed Time

Case 2C PP6523 Included Included Model 1 Case 2FCase 2BCase 2ACase

PP6523 Included Included Model 3 

PP702 Included Included Model 1 

PP702 Included Not included N/A 

In-depth Absorption Coefficient Model

In-depth Absorption

Radiative and Convective Heat Losses

PolymerResin

All cases are run with a external Incident heat flux of 30KW/m2

* Manually stopped because calculation become extremely slow when mass loss exceeds 50%
Impact of In-depth Absorption on Melting Proces

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-

Time interval between movie frames: 5 secondsThe inclusion of in-depth absorption leads to significant difference in shap

Melting process slows down Upper free surface moves much slower towards the back wallMelting layer becomes thicker

Temperature contour
Impact of Different Resin on Melting Process 

Case 2B (PP702, heat loss + in-depth absorption) Case 2C (PP6523, heat loss + in-d

Time interval between movie frames: 5 secondsPP6523 has a much lower viscosity at high temperature, which leads to faprocess

Temperature contour
Impact of Different Absorption Coefficient Model

Case 2B (PP702, in-depth absorption model 1)

Temperature contour

Case 2D (PP702, in-depth absorpti

Time interval between movie frames: 5 secondsModel 2 (Case 2D) has a much larger optical thickness than that in Model 

deeper penetration of incident heat flux below the surfacethicker melting layer and slower melting processlarge amount of material bulge over the bottom lip and causes the calculation to 
Impact of Different Absorption Coefficient Model

Case 2E (PP702, in-depth absorption model 3) Case 2F (PP6523, in-depth absorpt

Time interval between movie frames: 5 seconds

Temperature contour
Example: Velocity Vectors at t = 300 sec (zoom to 200% 

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-d

The addition of in-depth absorption leads to

slower melting process thicker melting layer below free surface

~50% mass loss ~15% mass loss
Mass Loss as Function of Time 

Case 2F 8.8* Computed from 250 s to 560 s

Case 2E 5.3* Computed from 260 s to 700 s

Case 2D 7.4* Computed from 340 s to 480 sCase 2C 9.0* Computed from 220 s to 440 s

0.15** Ohlemiller FAX 01-18-06 + Experiment

Case 2B 5.6* Computed from 220 s to 660 s

Case 2A 18.0* Computed from 200 s to 300 s, 

“steady state” Comments

Mass Loss Rate

Grid * p

** + m

2s

0

10203040506070

0 100 200 300 400 500 600 700 800

Time (s)

Mass Loss Percentage (%) Case 2A Case 2B Case 2C Case 2D Case 2E

Case 2F

Case 2A Case 2B

Case 2CCase 2D

Case 2E

Case 2F
0

0.020.040.060.08

0.1

0 0.01 0.02 0.0

Surface Velocity Magnitud

Elevation (m)

Case 2ACase 2BCase 2CCase 2ECase 2F

0

0.020.040.060.08

0.1

400 450 500 550 600 650 700 750

Surface Velocity Magnitude (m/s)

Elevation (m)

Case 2A, 280 secCase 2B, 440 secCase 2C, 350 secCase 2E, 485 secCase 2F, 380 sec

Example: Surface Temperature and Velocity Mag

35% mass l

35% mass loss

Surface temperature and velocity magnitude profiles at 35% mass 

Velocity magnitudes approach zero at upper no slip wallTemperature monotonically distributed along the surface in Case 2B, 2C, 2E, aInclusion of in-depth absorption decreases surface temperature and velocityOver-heating (T > 698K) occurs in Case 2A at upper portion of the free surface 

Measured surface temperature on a 2.5cmX25.0cm PP702 sample was 648 K (Ohlemiller FAX 01-18-06)
Summary

Studied Two Resins with Heat Losses and In-depth Absorption of Heat FluMost Cases were run to at least 50% mass loss

Calculation becomes extremely slow when mass loss exceeds 50%One case (Case 2D) diverged at 27% mass loss due to large optical thickness cOne case (Case 2E) is in progress using fitted absorption coefficient data

Computational Cost highlighted 

Less expensive than Task 1

Smaller computational domainLarger time step due to lower velocity caused by heat losses and in-depth absorption

Model features highlighted

Heat loss models and in-depth absorption modelsImpacts of different models on melting processesSolutions at selected timesVelocity field resolution across melt boundary layer

Example Results

Free surface location, temp., velocity - selected timesMovie of free surface location and temperatureMass loss, surface temperature, surface velocity, and comparison with experime
Investigating Impact of Different Moving 

Boundary Tracking Schemes on Predicte

Melting Behaviors

Contract No. SB1341-05-C-0041 

REACTION

INTERNATIONAL

ENGINEERING

06/22/2006

A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects
Outline

• Models for moving boundary tracking • Slumping problem with zero heat flux• Melting behaviors using different tracking models 
Figure 1

A

B

CD

F

No-slip adiabatic wall

Moving boundary

Ghost cell

E

Models for tracking moving boundary l

• The shape and location of movin

represented by grid point ACEF i

• The algorithm for determining mo

after each time step is described 

1. calculating new face center 

example

2. calculating new grid vertex 

example

3. Constrain grid vertex point 4. Smooth moving boundary: 

• To investigate the impact of the m

moving boundary location, two mcreated:

• Model 2: in step 

• Model 3: step 1-4 above t

constraint on wall point: 

• No-slip wall boundary affects the 

mass flux across the face center boundary, but has no direct impamovement of point A

A

m

x

x

CE

old

D

D

new

D

=

+

∆

ρ

&

x

y

( ) new

D

new

C

B

x

=

x

+

x

2

1

A

C

x

=

x

( 

A

C

new

x

x

x

C =

+

2

4

1

new

A B x

=

x

n

A

x
Slumping behavior as a function of vis

10

4 μ = 10 kg/m-s

3 kg/m-s 10

5 kg/m

VOF CPCFD VOF CPCFD VOF 

• External incident heat flux is set to zero in the CPCFD simulat• Model 1 for tracking moving boundary is used in all CPCFD sim
Vis = 1.e5, 

t = 555 s

Vis = 1.e4 

t = 55.5 s

Vis = 1.e3, 

t = 5.5 s

• External incident heat flux is set to zero in the CPCFD simulat• Model 1 of tracking moving boundary is used in all CPCFD sim• Similarity solutions have been achieved

Slumping behavior as a function of vis
Melting behavior with 30KW/m

2 Heat 

Model 1

Model 3

Model 2

t = 250 sec

Model 1

Model 2

Model 3

t = 350 sec t = 450 sec

Model 1

Model 2

Model 3

Model 1 

Model 3

PP6523

• The movement of the moving boundary at the top in Model 2 is significantly slowe

Model 1. This is because, in Model 2 (see figure 1), the movement of point A is appoint B rather than point C in model 1. Point B is halfway closer to the no-slip wa

• The different melting behaviors of Model 1 and Model 3 are caused by the surface 

process. 
Temperature Contour

T = 350 sec

Model 1 Model 2 Model 3

Contours drawn using only the interior cell center values 
Temperature Contour

T = 450 sec

Model 1 Model 2 Model 3

Contours drawn using only the interior cell center values 
Model 1 Model 3

PP6523
Melting behavior with 30KW/m

2 Heat 

0

10203040506070

0 200 400 600 800 10

Time (s)

Mass Loss Percentage (%) Model 1 Model 2

Model 3

PP6523

• Although the shape changes are very different (see the previous slide) when 

models to track the moving boundary, the differences in mass loss histories 

throughout the simulations.
Melting behavior with 30KW/m

2 Heat 

PP702N

• The different models for computing the location of the wall grid vertex point r

loss histories that are in close agreement for the duration of the simulations.

0

10203040506070

0 200 400 600 800 1000

Time (s)

Mass Loss Percentage (%) Case Model1 Case Model2

Case Model3
Parametric Study of Movin

Boundary Tracking Algorith

Qing Tang

REI

6/22/2006
Outline

• Under Relaxation (or Blending) of the Smoothing Function• Impact of Grid Refinement in the Vertical Direction• Impact of Under Relaxation of the Smoothing Function• Summary and Future Work* PP6523 is modeled in all cases
Figure 1

A

B

CD

F

No-slip adiabatic wall

Moving boundary

Ghost cell

E

x

y

• The algorithm for determining movin

each time step is described below

1. calculating new face center loc

example

2. calculating new grid vertex poi

example

3. Constrain grid vertex point to n4. Smooth moving boundary with 

5. re-apply constraint on wall poin

A

m

x x

CE

old D

D

new

D

= + ∆

ρ

&

( ) new D new C B x = x + x 21

A C x = x

( 

A 

old c new x x x 

C = ⋅ + − ⋅ 

41 α (1 α )

new new

A C

x = x

Under Relaxation of Smoothing Funct
Impact of Grid Refinement on Free Surf

300s 400s 500s 

Case 1: 30X80 GridCase 2: 30X160 Grida=0.0

Case 2 Case 2 Case 2

Case 1 Case 1 
Impact of Grid Refinement on Mass Loss 

0

10203040506070

0 100 200 300 400 500 600 7

Time (s)

Mass Loss Percentage (%) PP6523, 30X80 Grid

PP6523. 30X160 Grid

Case 1: 30X80 GridCase 2: 30X160 Grida=0.0
Under Relaxation Factor a=1.0, 0.99

t = 100 s

a=1.0a=0.99

t = 200 s t = 210 s t = 220 s 
Under Relaxation Factor a=0.9, 0.5, 0

t = 300 s t = 400 s t = 500 s t = 

a=0.9

a=0.5 a=0.0 a=0.0

a=0.9

a=0.5

a=0.9,0.5

a=0.0

a=0.9,0.5
Impact of Under Relaxation Factor on Mass L

0

10203040506070

0 100 200 300 400 500 600 

Time (s)

Mass Loss Percentage (%) alpha = 

alpha = alpha = 
Summary and Future Work

• Refining grid in the vertical direction does not have a significant im

shape of the moving boundary

• Refining grid in the vertical direction changes the mass loss rate s

the total mass loss reaches 30%, and makes the mass loss procesteady state 

• Smoothing the moving boundary is crucial to numerical stability.• Three different under relaxation factors in the smoothing function 

• Using small under relaxation factor, the moving boundary is 

different during the early melting stage.

• They all lead to similar moving boundary shapes when mass 

become significant

• smaller under relaxation factors lead to slightly faster mass lo

• Future Work: Test refined grid with small under relaxation factor
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

CFD Simulation of PP23K with steady hea

radiative and convective heat losses, an

depth absorption of radiative heat flu

(6/30/06)
PP23K, 30X80 Grid

0

10203040506070

0 30 60 90 120 150 180 210 240 270 300 

Time (s)

Mass Loss Percentage (%)

Mass Loss vs. Time
Temperature Contours

t = 100 sec. t = 300 sec.

T (K) 
Free Surface Evolution 

t = 25 sec t = 50 sec t = 100 sec t = 15

Velocity Magnitude (m/s)
Free Surface Evolution 

t = 200 sec t = 250 sec t = 300 sec t = 35

Velocity Magnitude (m/s)
Velocity Vectors Near Top Boundary

t = 25 sec

Velocity Magnitude (m/s)

t = 50 sec

Velocity Magnitude (m
Velocity Vectors Near Top Boundary

t = 75 sec

Velocity Magnitude (m/s)

t = 100 sec

Velocity Magnitude (m
Velocity Vectors Near Top Boundary

t = 125 sec

Velocity Magnitude (m/s)

t = 150 sec

Velocity Magnitude (m
Velocity Vectors Near Top Boundary

t = 200 sec

Velocity Magnitude (m/s)

t = 300 sec

Velocity Magnitude (m
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 2: 2D model with steady heat flux, ra

and convective heat losses, and in-dep

absorption of radiative heat flux

(7/11/06)
IntroductionRadiative and convective heat loss, In-depth absorption models

Only impacts energy equationBoth radiative and convective heat losses are only applied to the free surfacA heat source term is added to the energy equation to account for in-depth aAbsorption coefficient is a function of the distance from the free surface

Absorption coefficient models

In-depth absorption coefficient data for one polymer resin provided by NIST

Email sent by K. Butler to M. Bockelie on 2/22/2006Data applicable to PP23K in the HIFT apparatus

Cases studied

Computational gridCase summaryExample results: 

Free surface location, temp., velocityMass loss, surface temperature, surface velocity, and comparisons with data

Summary

Outline
Several modeling issues have been investigated during the perfoperiod of Task 2, they are

In-depth absorption modelImplementation of top no-slip wall boundaryAdiabatic wall boundaryGrid smoothing functionGrid sensitivity

Simulation of a third polymer raisin (PP23K) has also been perfoResults and discussions about these investigations can be found PowerPoint progress reports and e-mails with NIST during May 2006Results shown in this report represent the latest from simulating raisin PP702 and PP6523 based on the experience gained throu

Introduction
Heat loss and In-depth absorption models

I0

qrad

I

Figure 1, Sketch of in-depth 

absorption and heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy equation source terms due to in-depth absorption, radiative and convect

In-depth absorption: , assuming absorption process in x-direction oRadiative heat loss: , assuming radiative heat loss only on fConvective heat loss:Parameters: A – surface area; - volume; - emissivity (=1.0); - Stefan-Boltzmann Consta

T0 – ambient temperature (=298 K); hconv – convective heat transfer coefficient (= 

A heat sink term is added to the energy equation to prevent overheating (i.e., T > 750K) The models are implemented in subroutine fblck_bc_sca

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

σ
Absorption Coefficient Model

Figure 2, Incident heat flux distribution in the horizontal direction 

In-depth absorption model 

⎪⎩⎪⎨⎧

=

≤

<

=

≤

<

− − +− +

I/I 10 , 0.3

x 1.3 (cm)

I/I 10 , 0

x 0.3 (cm)

( 0.5154 1.473x 0.4339x )

0

( 5.17x 6.964x )

0

2

2

0

0.20.40.60.8

1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x, Distance from Free Surface (cm)

Normalized Incident Heat Flux I/I 0
Computational Grid

Polymer sample dimension: 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Figure 3, Computational grid at t = 0 sec. 
Case Summary

0.90.00.0Surface Smoothing Underrelaxation Factor ( )*

~ 28 h~ 26 h~ 50 hEstimCPU Time*

70%, 700 sec65%, 700 sec60%, 370 secMass Loss and Simulated Time

Case 2C PP6523 Included IncludedCase 2BCase 2ACase

PP702 Included Included

PP702 Included Not included

In-depth Absorption

Radiative and Convective Heat Losses

PolymerResin

External Incident heat flux: 30KW/m2

* As per PowerPoint progress report 06-22-2006: **all cases were run on a 3.0GHz Intel Xeon workstation

( 

old A old C new Cx = ⋅ x + − ⋅ x 

41 α (1 α )

α
Impact of In-depth Absorption on Melting Proces

Case 2A 

(PP702, heat loss only)

Case 2B 

(PP702, 

heat loss + 

in-depth absorption)

Time interval between movie frames: 5 secondsThe inclusion of in-depth absorption leads to significant difference in shape chan

Melting process slows down Upper free surface moves much slower towards the back wallMelting layer becomes thicker

Temperature contour

Case 2C 

(PP6523, 

heat loss + 

in-depth absorpti
Example: Temperature Contours at t = 250 and 65

Case 2B (PP702, heat loss + in-depth absorption)Contour lines are perpendicular to the solid walls in the near wall regions, which indicates the walls are adiabatic 

t = 250 sec t = 650 sec
Example: Velocity Vectors at t = 300 sec

Case 2A (PP702, heat loss only) Case 2B (PP702, heat loss + in-d

The addition of in-depth absorption leads to

slower melting process thicker melting layer below free surface

~51% mass loss ~19% mass loss

Vectors in different boxes are not the same scale
Mass Loss as Function of Time 

Case 2C 7.5* Computed from 240 s to 600 s

Experiment 0.15** Ohlemiller FAX 01-18-06

Case 2B 6.5* Computed from 200 s to 600 sCase 2A 24.0* Computed from 200 s to 300 s, 

“steady state” Comments

Mass Loss Rate

Grid

*

0

10203040506070

0 100 200 300 400 500 600 700

Time (s)

Mass Loss Percentage (%) Case 2A Case 2B

Case 2C
0

0.020.040.060.08

0.1

0 0.01 0.02 0.03 0.04 0.05 0.06 

Surface Velocity Magnitude (

Elevation (m)

Case 2A, 255 secCase 2B, 382 secCase 2C, 365 sec

0

0.020.040.060.08

0.1

500 525 550 575 600 625 650 675 700 725 750

Surface Temperature (K)

Elevation (m)

Case 2A, 255 secCase 2B, 382 secCase 2C, 365 sec

Example: Surface Temperature and Velocity Mag

35% mass loss 35% mass loss

Surface temperature and velocity magnitude profiles at 35% mass lo

Velocity magnitudes approach zero at upper no slip wallInclusion of in-depth absorption decreases surface temperature and velocity mOver-heating (T > 698K) occurs in Case 2A where surface temperature is at the near the top 

Measured surface temperature on a 2.5cmX25.0cm PP702 sample was 648 K (Ohlemiller FAX 01-18-06)
Summary

Results presented for Two Resins with Heat Losses and In-depth AbsorptiCases were run to at least 60% mass lossComputational Cost highlighted 

Less expensive than Task 1

Smaller computational domainLarger time step due to lower velocity caused by heat losses and in-depth absorption

Model features highlighted

Heat loss model and in-depth absorption modelImpacts of different models on melting processesSolutions at selected timesVelocity field resolution across melt boundary layer

Example Results

Free surface location, temp., velocity - selected timesMovie of free surface location and temperatureMass loss, surface temperature, surface velocity, and comparison with experime

CPCFD source code for Task 2 has been provided to NIST via Email (6/30
Task 3
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 3: 2D model with in-depth gasifica

(8/11/06)
IntroductionIn-depth gasification modelsCases studied

Computational gridCase summaryExample results: 

Free surface location, tempMass loss, surface temperature, and comparisons with data

Summary

Outline
In-depth gasification model has been implemented into the CPCGasification is represented by the removal of polymer mass anusing sink terms in both mass and energy equations.Simulations of two polymer raisins (PP702N and PP6523) have performedImpact of in-depth absorption on gasification has been investigResults shown in this report represent the latest from simulating raisin PP702N and PP6523 based on the experience gained thand Task 2.

Introduction
Gasification, absorption and heat losses models

I0

qrad

I

Figure 1, Sketch of gasification, 

absorption and heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy source/sink terms due to gasification, absorption, and heat losses

In-depth gasification:In-depth absorption: , assuming absorption process in x-direction oRadiative heat loss: , assuming radiative heat loss only on fConvective heat loss:Heat sink term to avoid over heating in Task 1 and 2 has been deactivated in this study

Mass sink terms due to gasification: 

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

σ

( )

∫

=

−

Ω

-E/RT

q& gasification

HvρAe

d

Ω

( )

∫

=

−

Ω

-E/RT

m& gasification

ρAe

d

Ω

qgas, mgas
Gasification Model Parameters

Figure 2, Mass loss rate due to gasification as functions of time 

05

101520

350 450 550 650 750

Temperature (K)

Gasification rate (kg/s m3)

PP6523PP702N

PP6523 2.23e+13 26000 8.0e+5 900PP702N 2.18e+12 24400 8.0e+5 900

r

(kg/m3

Hv

(J/kg)

A E/R

(1/s)

Polymer Raisin

Gasification model parameters
Computational Grid

Polymer sample dimension: 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Figure 3, Computational grid at t = 0 sec. 
Case Summary

includedIncludedIncludedIncludedIn-depth Gasification

Case 3C PP6523 Included Included 65%, 600 sec ~ 24 h

~ 30 h~ 36 h~ 24 hEstimateCPU Tim

60%, 280 sec55%, 230 sec65%, 600 secMass Loss and Simulated Time

Case 3D PP6523 Included Not includedCase 3BCase 3ACase

PP702N Included Not included

PP702N Included Included

In-depth Absorption

Radiative and Convective Heat Losses

PolymerResin

External Incident heat flux: 30KW/m2
Simulated Melting Process (PP702N)

Case 3A 

(PP702N, with in-depth absorption)

Time interval between movie frames: 5 seconds

Temperature 

contour

Case 3B 

(PP702N, without in-depth absorp
Simulated Melting Process (PP6523)

Case 3C 

(PP6523, with in-depth absorption)

Time interval between movie frames: 5 seconds

Temperature 

contour

Case 3D 

(PP6523, without in-depth absorpt
Gasification Rate (kg/s) (PP6523)

Case 3C 

(PP6523, with in-depth absorption)

Case 3D, 255 sec. 

(PP6523, without in-depth absorp

T (K) 

T (K)

GasificationRate (kg/s)
0

10203040506070

150 170 190 210 230 250 270

Time (s)

Mass Percentage (%)

Case 3B: Total Mass LossCase 3B: Mass Gained by Catch Pan

0

10203040506070

0 100 200 300 400 500 600 700

Time (s)

Mass Percentage (%)

Case 3A: Total Mass LossCase 3A: Mass Gained by Catch Pan

Mass Loss as Function of Time (PP702N) 

• The differ

between tand the grepresent loss due t

• Mass loss 

gasificationegligible In Case 3mass loss gasificatioof the sim

• Experime

indicates gasificatiofor 8 to 10

Case 3B

Case 3A
Mass Loss as Function of Time (PP6523) 

0

10203040506070

0 100 200 300 400 500 600 700

Time (s)

Mass Percentage (%)

Case 3C: Total Mass LossCase 3C: Mass Gained by Catch Pan

0

1020304050607080

150 170 190 210 230 250 270 290

Time (s)

Mass Percentage (%)

Case 3D: Total Mass LossCase 3D: Mass Gained by Catch Pan

Case 3D

Case 3C

• The diff

between and the represeloss due 

• Mass lo

gasificanegligiband is vCase 3D end of t
0

0.020.040.060.08

0.1

1.00E-09 1.00E-08 1.00E-07 1.00E-06 1

Magnitude of gasification sink term near surf

Elevation (m)

Case 3ACase 3B

0

0.020.040.060.08

0.1

300 350 400 450 500 550 600 650 700 750 800

Surface Temperature (K)

Elevation (m)

Case 3ACase 3B

Example: Surface Temperature and Gasification 

35% mass loss 35% mass loss

Surface temperature and gasification sink term profiles at 35% mass loss 

Surface temperature increases more then 100 K when in-depth absorption moMagnitude of gasification sink term near the surface decreases by three order

when in-depth absorption model turned off

Predicted surface temperature is in good agreement with measurement when absorption model turned on

Measured surface temperature on a 2.5cmX25.0cm PP702 sample was 648 K (Ohlemiller FAX 01-18-06)
0

0.020.040.060.08

0.1

1.00E-09 1.00E-08 1.00E-07 1.00E-06 1.

Magnitude of gasification sink term near surfa

Elevation (m)

Case 3CCase 3D

0

0.020.040.060.08

0.1

300 350 400 450 500 550 600 650 700 750 800

Surface Temperature (K)

Elevation (m)

Case 3CCase 3D

Example: Surface Temperature and Gasification 

35% mass loss 35% mass loss

Surface temperature and gasification sink term profiles at 35% mass loss 
Summary

Results presented for two raisins with gasification, heat losses and in-deptexternal heat flux at 30 KW/m2Cases were run to at least 60% mass lossComputational Cost highlighted 

Inclusion of gasification model does not have a big impact on run time.

• Example Results

Free surface location, temp., gasification sink termMovie of free surface location and temperatureMass loss, surface temperature, gasification sink term 

Discrepancy with measurement

Measured surface temperature and gasification mass loss for PP702N are availaPredicted surface temperature agrees well with measurement when all sub-mod

gasification, heat loss, and in-depth absorption models, are included in the simu

predicted gasification mass loss is negligible.

When turning off the in-depth absorption model, predicted gasification mass loss significant. However, the predicted surface temperature is at least 100 K higher t

measurement 

The above observations suggest that, in order to match the experimental data, b

absorption model and the in-depth gasification model need improvments. 
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 3: 2D model with in-depth gasifica

Parametric Study of model constant

(8/31/06)
IntroductionIn-depth absorption and gasification models in the parametric stParametric Cases studied

Increase gasification rateDecrease heat flux penetration lengthShift gasification to start at a lower temperatureIncrease external heat flux

Summary

Outline
In-depth absorption and gasification models provided by NIST himplemented into the CPCFD code. See 8/11/06 report for detaAs noted in 8-11-06 report, after simulating two polymer resins PP6523), it was found that

Predicted surface temperature agrees well with measurement when all sub-m

gasification, heat loss, and in-depth absorption models, are included in the siHowever, the predicted gasification mass loss is negligible.

When turning off the in-depth absorption model, predicted gasification mass significant. However, the predicted surface temperature is at least 100 K high

measurement

The above observations suggest that, in order to match the expdata, both the in-depth absorption model and the in-depth gasifneed adjustments.A parametric study regarding the two models is summarized in Simulations are performed for resin PP702N only

Introduction
Gasification, absorption and heat losses models

I0

qrad

I

Sketch of gasification, absorption and 

heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy source/sink terms due to gasification, absorption, and heat losses

In-depth gasification:In-depth absorption: , assuming absorption process in x-direction Radiative heat loss: , assuming radiative heat loss only on fConvective heat loss:Heat sink term to avoid over heating in Task 1 and 2 has been deactivated in this study

Mass sink terms due to gasification: 

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

( )

∫

=

−

Ω

-E/RT

q& gasification

HvρAe

d

Ω

( )

∫

=

−

Ω

-E/RT

m& gasification

ρAe

d

Ω

qgas, mgas
Gasification Models Used in Parametric Study

Model G2 2.18e+13 Model G3 2.18e+12 

Model G4 2.18e+12 Model G5Model G1

2.18e+12 2.18e+12 

A

(1/s)

Gasification model 

Mass loss rate due to gasification as functions of 

time 

=

Model expression:

y

x: temperature (K)y: gasification rate (kg/s 

Specific gasification heat increases from 800 J/g to 1260 J/g for calculating the gasif

term in all models

Model G1 is the base modelModel G2 increases the magnitude of gasification rate by a factor of 10Model G3, G4 and G5 shift the occurrence of gasification to lower temperature by 10

respectively.

0

500

100015002000250030003500

300 400 500 600 700 800 900

Temperature (K)

Gasification Rate (kg/s m3)

Model G1Model G2Model G3Model G4Model G5
Absorption Models Used in Parametric Study

Normalized Incident Heat Flux Profile as a 

Function of Distance from Free Surface

Model expression:

Model A1:Model A2:

Model A1 is the base modelModel A2 shortens the penetration length of the incident heat fluxThe heat flux penetration length covers at least 4 cells in the horizontal direction nea

surface

( )( 

10 10 , 0

2

0.5154 1.473 0.4339

0

5.17 6.964

0

==

− − +− +

I

I

I

I

x 

x x

0

0.20.40.60.8

1

1.2

0 0.2 0.4 0.6 0.8 1 1.2

Distance (cm)

I/I0

Model A1Model A2

exp( 153.507

7.6

2

0

=

−

x

−

I

I

x: distance from free surface in 
Computational Grid

Polymer sample dimension: 5.0cm X 10.0cm X 2.54cmModeled as a 2-D problem with a computational domain of 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Computational grid at t = 0 sec. 
Case Summary

Case 3K Model G3 Model A2 40Case 3J Model G5 Model A2 30

Case 3L Model G5 Model A2 40

Model A2Model A2Model A2Model A1Model A1

In-depth Absorption Model

3030303030

External heat flux(KW/m2)

Case 3I Model G4Case 3H Model G3

Model G1Model G2Model G1

In-depth Gasification Model

Case 3GCase 3FCase 3ACase #

All cases are for PP702NRadiation and convective heat losses at free surface are included in all Case 3A is the baseline case
Simulation Results Summary

Case 3K 602 640 2.93e-5 0.080 

648 N/A 0.053 

Measurement N/A ! * +

Case 3L 355 654 1.92e-5 0.15 

Measurement * N/A N/A N/A 0.086 

+

Case 3J 675 637 9.10e-6 0.083 

Case 3I 461 646 4.15e-6 0.15 Case 3H 806 625 1.53e-5 0.063 

Case 3G 378 649 3.47e-7 0.18 

Case 3F 383 634 1.44e-6 0.17 

Case 3A 368 636 1.54e-7 0.18 

M

Total mass 

loss rate*

(g/s)

Average surface 

gasification rate at 

35% mass loss 

(g/s m)

Average surface 

temperature at 

35% mass

loss (K)

Exposed time to 

reach 35% 

mass loss (s)

* mass loss rates are calculated from 200s to 600s for cases with 30KW/m2 heat flux and from 100 to 

with 40KW/m2 heat flux. The initial dimension of the simulated sample is 5.0cmX10.0cmX2.54cm. + measurement data was calculated based on data reported in “A Progress Report on Numerical Mod

Experimental Polymer Melt Flow Behavior” by Kathryn M. Butler, et al. The initial dimension of the 

sample was not reported in the report and is assumed to be the same as that of the computa

! Measured surface temperature was for a 2.5cmX25.0cm PP702 sample (Ohlemiller FAX 01-18-06). 

was a point measurement (i.e., not a measurement of the average surface temperature). At what timtime period that this measurement was done was not available.
Impact of Increasing Gasification Rate

Mass gained by catch pan as a function of time

Gasification mass loss as a function of time

0

1020304050607080

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3F

0

0.10.20.30.40.50.60.70.80.9

1

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3F

0

0.020.040.060.08

0.1

600 610 620 630 640 650 6

Surface Temperature (K)

Height (m)

Case 3ACase 3F

Surface Temperature Profile at 35% m

Increasing gasification rate by a fadecreases surface temperature slGasification mass loss in Case 3F that of the baseline caseCompared to mass loss due to megasification mass loss is still negliOverall mass loss rate in Case 3F smaller than that of Case 3A due tlower surface temperature and thastill dominates the mass loss.
Impact of Decreasing Heat Flux Penetration Leng

Mass gained by catch pan as a function of time

Gasification mass loss as a function of time

0

0.05

0.1

0.15

0.2

0.25

0.3

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3G

0

1020304050607080

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3G

Surface Temperature Profile at 35% m

0

0.020.040.060.08

0.1

600 610 620 630 640 650 660

Surface Temperature (K)

Height (m)

Case 3ACase 3G

Decreasing incident heat flux peneincreases surface temperatureIncreased surface temperature leagasification mass lossCompared to mass loss due to megasification mass loss is still negliDecreasing heat flux penetration ldoes not have a big impact on ovehistory
Impact of Shifting Gasification to Lower Tempera

Mass gained by catch pan as a function of time

Gasification mass loss as a function of time

0

1020304050607080

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3HCase 3ICase 3J

02468

10

0 200 400 600 800

Time Exposed (s)

Mass (g)

Case 3ACase 3HCase 3ICase 3J

0

0.020.040.060.08

0.1

600 610 620 630 640 650 6

Surface Temperature (K)

Height (m) Case 3A

Case 3HCase 3ICase 3J

Surface Temperature Profile at 35% m

In Case H, I and J the modified abis adaptedShifting gasification to lower temp

has a big impact on the melting pr

Larger temperature shift (value of 

leads to lower surface temperatureloss due to melting flow, but more to gasification
Close Look at Case H, I and J

620625630635640645650

25 50 75 100 125

Temperature shift (

C ) in gasification model (K)

Surface teperature at 35% mass loss (K) Simulation Measurement

0.0E+00

4.0E-068.0E-061.2E-051.6E-052.0E-05

25 50 75 100 125

Temperature shift (

C ) in gasification model (K)

Gasification rate at 35% mass loss (g/s m)Simulation 

0

0.020.040.060.08

0.1

0.120.140.16

25 50 75 1

Temperature shift (

C ) in gasificatio

Surface teperature at 35% mass loss (K)

Simulation: total mass lossMeasurement: total mass lossSimulation: mass loss due to meMeasurement: mass loss due to Simulation: mass loss due to gasMeasurement: mass loss due to 

Case 3H (C=100K) seems to yield agreement with the measured data mass loss rates
Mass Loss Histories at Different Heat Flux

Symbols – experimental data (“A Progress Report on Numerical Modeling of Experim

Flow Behavior” by Kathryn M. Butler, et al)

Lines – Case 3H and Case 3K Simulation results

0

102030405060

0 200 400 600 800

Time Exposed (s)

Mass (g)

30KW, total mass loss30KW, melting flow mass loss 40KW, total mass loss40KW, melting flow mass loss 30KW, total mass loss30KW, melting flow mass loss 40KW, total mass loss40KW, melting flow mass loss
Summary

The impacts of different in-depth absorption and gasification models on the 

behavior have been studied

When only modifying one model, we have the following observation:

Increasing the magnitude of the gasification rate by a factor of ten significantly e

gasification mass loss, which, however, is still negligible to the mass loss due to 

Shortening the penetration length of incident heat flux in the in-depth absorption 

of at least three increases the surface temperature by about 15K, which leads to 

gasification mass loss. However, the overall mass loss history does not change 

When modifying both models, we found that:

Shifting the occurrence of gasification towards lower temperature range significa

melting behavior.

More temperature shifting towards lower temperature range -> lower surface tem

mass loss due to melting flow -> more mass loss due to gasification -> slower ov

process.

When combining model G3 for gasification and model A2 for absorption, th

results for PP702N seem to agree well with experimental data for heat flux 

40KW/m2.
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 3: 2D model with in-depth gasifica

Parametric Study of model constant

(9/14/06)
Physical ModelCFD ModelComparison of previous simulation results with new experimentIn-depth gasification models to be tested

Outline
Physical Model

10 cm 5 cm

10 cm

Radiation

Feed Direction

Melting Front
CFD Model

5 cm

10 cm

Radiation

No-slip walls

-u
Comparison of CFD Results with Data

0.050.100.150.200.250.300.35

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Total Mass Loss Rate (g/s) Measurement

Simulation

Simulation results are from Case 3H and Case 3K using CFD model 1 (see 8/31/06 rMeasurement data provided by NIST (TOM, O., Fax 9/8/06)
Comparison of CFD Results with Data

Simulation results are from Case 3H and Case 3K using CFD model 1 (see 8/31/06 rMeasurement data provided by NIST (TOM, O., Fax 9/8/06)

0.00.20.40.60.81.0

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Fraction of Mass Loss Occuring as Melt Flow Measurement Simulation
Comparison of CFD Results with Data

Simulation results are from Case 3H and Case 3K using CFD model 1 (see 8/31/06 rMeasurement data provided by NIST (TOM, O., Fax 9/8/06)

300350400450500

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Surface Temperature (°C)

MeasurementSimulation
Gasification Models for Future Work

Model G2 2.18e+13 

Model G6 3.27e+14 Model G1 

2.18e+12

A

(1/s)

Gasification model p

Mass loss rate (in log scale) due to gasification as 

functions of time 

⎟⎠⎞ ⎜

⎝⎛

=

−

x

B Model expression:

y

Aexp

x: temperature (K)y: gasification rate (kg/s m

3) 

Model G1 is the base modelModel G2 increases the magnitude of gasification rate by a factor of 10Model G6 increases the magnitude of gasification rate by a factor of 150

1.E-211.E-171.E-13

1.E-91.E-51.E-1

1.E+31.E+7

300 400 500 600 700 800 900

Temperature (K)

Gasification Rate (kg/s m3)

Model G1Model G2Model G6

0.E+05.E+31.E+42.E+42.E+4

300 400 500 600 700 

Temperature (K)

Gasification Rate (kg/s m3)

Model G1Model G2Model G6

Mass loss rate due to gasification as f

time 
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 3: 2D model with in-depth gasifica

Parametric Study of model constant

(9/26/06)
Physical ModelCFD ModelModifications in absorption and gasification modelsSelected simulation results and comparisons with data for PP70Selected simulation results for PP6523Summary

Outline
Physical Model

W = 10 cm D = 5 cm

H = 10 cm

Radiation

Feed Direction

Melting Front

2D “slice” for CFD model

h x w x d = 10cm x 10cm x 5 cm• in CFD model w only used to compare to test da
Computational Grid

• Modeled as a 2-D problem 

with a computational domain of 5.0cm X 10.0cm

• 30X80 Grid (biased)• Initial grid resolution the same 

as the 19X200 grid used in Task 1

• Equivalent to 100X80 grid if 

using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Computational grid at t = 0 sec. 

D = 5 cm

H = 10 cmRadia
Gasification Models

Model G6 3.27Model G1 

2.18

A 

(1

Gasification mode

Mass loss rate (in log scale) due to gasification as 

functions of time 

⎟⎠⎞ ⎜⎝⎛ = − xB y Aexp

Model express

x: temperature y: gasification 

Model G1 is the base modelModel G6 increases the magnitude of gasification rate by a factor of 150

0.0E+05.0E+31.0E+41.5E+42.0E+4

300 400 500 600 700 800

Temperature (K)

Gasification Rate (kg/s m

3

)

Model G1Model G6
Absorption Models

Normalized Incident Heat Flux Profile as a 

Function of Distance from Free Surface

Model expression:

Model A1:Model A3:

Model A1 is the base modelModel A3 shortens the penetration length of the incident heat fluxThe heat flux penetration length covers at least 4 cells in the horizontal direction nea

surface

( )( 

10 

10 , 0 

2

0.5154 1.473 0.4339

0

5.17 6.964

0

= 

= 

− − +− +

I

I

I

I

x 

x x

exp( 11.513 ), 

0

= − x 

I

I

x: distance from free surface in 

0

0.20.40.60.8

1

1.2

0 0.2 0.4 0.6 0.8 1 1.2

Distance (cm)

I/I

0

Model A1Model A3
Case Summary (PP702N)

Model A3Model A3Model A3Model A1

In-depth Absorption Model

20463030

External heat flux(KW/m2)

Case 3Q Model G6

Model G6Model G6Model G1

In-depth Gasification Model

Case 3PCase 3NCase 3ACase #

All cases are for PP702NRadiation and convective heat losses at free surface are included in Case 3A is the baseline caseModel G6 and A3 in Case 3N are developed to match measured masat 30kW/m

2 heat flux

Case 3P and Case 3Q use same model parameters as Case 3N
0

50

100150200250300

0 200 400 600 800

Time Exposed (s)

Mass Loss (g)

Impacts of Changing Gasification and Absorption 

Case 3A (fraction of mass as melt flow > 99

Case 3N

Mass loss history 

- solid lines = total mass loss- dashed lines = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasifica
Impacts of External Heat Flux (PP702N)

0

50

100150200250

0 200 400 600 800

Time (s)

Mass loss (g)

Case 3P

Case 3N

Case 3Q

Mass loss history 

- solid lines = total mass loss- dashed lines = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasifica
Simulation vs. Measurement (PP702N)

0.000.150.300.450.60

10 20 30 40 

Incident Radiant Flux (kW/m2)

Total Mass Loss Rate (g/s) Measurement Simulation, Case 3N, 3P & 3Q Simulation, Case 3A
Simulation vs. Measurement (PP702N)

0.00.20.40.60.81.0

10 20 30 40 

Incident Radiant Flux (kW/m2)

Fraction of Mass Loss Occuring as Melt Flow Measurement Simulation, Case 3N, 3P & 3

Simulation, Case 3A
Simulation vs. Measurement (PP702N)

0

100200300400500

10 20 30 40 

Incident Radiant Flux (kW/m2)

Surface Temperature (°C) Measurement

Simulation, Case 3N, 3P & Simulation, Case 3A
Gasification Models for PP6523

0.0E+05.0E+31.0E+41.5E+42.0E+42.5E+4

300 400 500 600 700 800

Temperature (K)

Gasification Rate (kg/s m3)

Model G7Model G8

Model G8 3.Model G7 

2.

Gasification mo

⎟⎠⎞ ⎜

⎝⎛

=

−

x

B

y

Aexp

Model expres

x: temperatuy: gasificatio

Mass loss rate (in log scale) due to gasification as 

functions of time 

Model G7 is the base modelModel G8 increases the magnitude of gasification rate by a factor of 150Parameter for Model G8 is selected based on the experience of modeling PP702N
Case Summary (PP6523)

3030

External heat flux(KW/m2)

Model A3Model A1

In-depth Absorption Model

Model G8Model G7

In-depth Gasification Model

Case 3ICase 3CCase #

Radiation and convective heat losses at free surface are included in Case 3C is the baseline case
0

50

100150200250300

0 500 1000 1500 2000

Time Exposed (s)

Mass Loss (g)

Predicted Mass Loss History for PP6523

Case 3C (fraction of mass loss occurring as melt flow > 99%)

Case 3I

Mass loss history 

- solid lines = total mass loss- dashed lines = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasifica
Case Summary (PP6523)

3030

External heat flux(KW/m2)

0.065

0.69

Total mass loss rate (g/s)

0.280.99

Fraction of mass loss due to melt flow

Model A3Model A1

In-depth Absorption Model

SurfaTemp25% 

Model G8Model G7

In-depth Gasification Model

Case 3ICase 3CCase #

Radiation and convective heat losses at free surface are included in Case 3C is the baseline case
Summary

The in-depth absorption and gasification models are modified to match measured total mass loss rate and the measured fraction of mass loss as melt flow at 30 kW/m2 external heat flux for PP702NCase 3A (baseline) uses the absorption and gasification models provid

NIST does not match measured data for mass loss or surface tempera

Case 3N which uses the modified models shows 

3⁄4 vastly different simulation results (compare to baseline)

3⁄4 matches the measured data reasonably well.

The same modified models developed for PP702N exposed to 30kW/m

have been applied to cases with 20 and 46 kW/m2 heat flux. The predic

results agree reasonably well with measured data and show correct tretotal mass loss rate, fraction of mass loss due to melt flow, and surface temperature as a function of heat flux.

Two simulations have also been performed for PP6523 exposed to 30 heat flux using different in-depth absorption and gasification models. It that the results are highly sensitive to the model parameters in the mod
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 3: 2D model with in-depth gasifica

(12/05/06)
IntroductionIn-depth gasification modelsIn-depth absorption modelsCases studied

Computational gridCase summaryExample results: 

Free surface location, temperature contourMass loss, surface temperature, and comparisons with data

Summary

Outline
In-depth gasification model has been implemented into the CPCGasification is represented by the removal of polymer mass anusing sink terms in both mass and energy equations.Simulation results of PP702N are presented in this reportA new in-depth absorption model is used in the simulations.

Introduction
Gasification, absorption and heat losses models

I0

qrad

I

Figure 1, Sketch of gasification, 

absorption and heat loss model

no-slip wall

Polymer Sample(5.0cmX10.0cm)

Free surface

qconv

y

x

Energy source/sink terms due to gasification, absorption, and heat losses

In-depth gasification:In-depth absorption: , assuming absorption process in x-direction onRadiative heat loss: , assuming radiative heat loss only on freConvective heat loss:Heat sink term to avoid over heating in Task 1 and 2 has been deactivated in this study

Mass sink terms due to gasification: 

∫

=

Ω absorption

d

Ω

dx

dI

q&

q

A

σ(T

T )

4

0

4 &

radiation

=

−

ε

−

q Ah (T

T ) convection

=

−

conv

− 0 &

Ω

ε

σ

( )

∫

=

−

Ω

-E/RT

q& gasification

HvρAe

d

Ω

( )

∫

=

−

Ω

-E/RT

m& gasification

ρAe

d

Ω

qgas, mgas • ρ = 90

• cp = 24

• k = 0.2• q0 = 30 

• T0 = 29

• ε = 1.0• σ = 5.6

• hconv = 

• Hv = 1
Gasification Model Parameters

Gasification mode

⎟⎠⎞ ⎜

⎝⎛

=

−

x

B Model expression:

y

ρAexp

x: temperature (K)y: gasification rate (kg/s m

3) 

Model G1 is the base modelModel G6 increases the magnitude of gasification rate by a factor of 150G1 and G6 are used in the current calculations

Mass loss rate (in log scale) d

function of tim

1.E-211.E-171.E-13

1.E-91.E-51.E-1

1.E+31.E+7

300 400 500 600 

Temperature (

Gasification Rate (kg/s m3)

0.E+05.E+31.E+42.E+42.E+4

300 400 500 600 700 800 900

Temperature (K)

Gasification Rate (kg/s m3)

Model G1Model G2Model G6

Mass loss rate due to gasification as function of time 

Model G2 2.18e

Model G6 3.27eModel G1 

2.18e

A

(1/s
In-Depth Absorption Model

from Tom Ohlemiller 10/25/06

Model A4:
Physical Model

W = 10 cm D = 5 cm

H = 10 cm

Radiation

Feed Direction

Melting Front

2D “slice” for CFD model

h x w x d = 10cm x 10cm x 5 cm• in CFD model w only used to compare to test da
Computational Grid

Modeled polymer sample dimension: 5.0cm X 10.0cm30X80 Grid (biased)Initial grid resolution the same as the 19X200 grid used in Task 1Equivalent to 100X80 grid if using uniform grid spacing and fixed grid

5.e

No-slip wall

Free surf

Figure 3, Computational grid at t = 0 sec. 
Case Summary

G1G6In-depth Gasification

~ 10 h~ 20 hEstimateCPU Tim

65%, 600 sec55%, 1000 secMass Loss and Simulated Time

Case 3SCase 3RCase

PP702N Included A4

PP702N Included A4

In-depth Absorption

Radiative and Convective Heat Losses

PolymerResin

External Incident heat flux: 30KW/m2
Simulated Melting Process (PP702N)

Case 3R 

(0 ~ 1000 s, A4 + G6)

Time interval between movie frames: 5 seconds

Case 3S 

(0 ~ 600 s, A4 + G1)

Temperature 

contour
Mass Loss as Function of Time

Mass loss history 

- solid lines = total mass loss- dashed lines = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasifica

0

10203040506070

0 200 400 600 800

Time (s)

Mass Loss Percentage (%)

Case 3RCase 3RCase 3SCase 3S

0

50

100150200250300350

0 200 400 600 

Time (s)

Mass Loss (g)

Case 3RCase 3RCase 3SCase 3S
Simulation vs. Measurement

0.000.150.300.450.600.75

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Total Mass Loss Rate (g/s) Measurement

Simulation, Case 3RSimulation, Case 3S

Case 3R, rate (0.28 g/s) computed using data from 110 s to Case 3S, rate (0.74 g/s) computed using data from 210 s to 
Simulation vs. Measurement (PP702N)

0.00.20.40.60.81.0

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Fraction of Mass Loss Occuring as Melt Flow Measurement Simulation, Case 3R Simulation, Case 3S

Case 3R, fraction (0.997) computed using data of 900 sCase 3S, fraction (0.740) computed using data of 600 s 
Simulation vs. Measurement (PP702N)

250300350400450500

10 20 30 40 50

Incident Radiant Flux (kW/m2)

Surface Temperature (°C) Measurement

Simulation, Case 3RSimulation, Case 3S

Predicted surface temperatures are area-averaged temperaturefree surface at 450 s and 300 s for Case 3R (358°C) and Case respectively.
Summary

Results presented for PP702N with gasification, heat losses and in-depth aexternal heat flux at 30 KW/m2Cases were run to at least 55% mass lossComputational Cost highlightedExample Results

Movie of free surface location and temperatureMass loss rate, fraction of mass loss occurring as melt flow and surface tempera

• Observations

Free surface near the top boundary recedes faster initiallyFree surface curvature near the bottom lip is suppressed artificially by the smoot

to stabilize the computation

Increasing the gasification rate by a factor of 150 from model G1 to G6 leads to 

Lower surface temperature (i.e., decreasing by about 20°C)Lower total mass loss rateHigher fraction of mass loss occurring as melting flow

Comparison with data

Using the original gasification model G1, the predicted total mass loss rate is fou

the measurement; the predicted fraction of mass loss due to gasification is negli

smaller than the measurement; though predicted surface temperature appears in 

with the data.

Using the modified gasification model G6, the predicted total mass loss rate is 5

measurement; the predicted fraction of mass loss occurring as melt flow and sur

are in good agreement with the data.
Task 4
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

Goals, Methodology, and Plan

(01/10/07)
Task descriptionModeling methodologyPreliminary plan

Outline
Add to the existing model the flow of melt on the solid catch su

Include catch surface material properties (slip/non-slip, surface temFinite thickness of catch basinEffects of surface tension Nonzero gap between sample and catch surfaceNonzero tilt angle for the catch surfaceNo change to in-depth absorption and in-depth gasification model

Task 4A – model PP702 with specified heat flux, gap and tilt anTask 4B – model PP6523 (same conditions as those used in TaTask 4C – model PP702 or PP6523 with a different gap (no chaconditions)Task 4D – NIST site visit (project update)Task 4E – NIST site visit (base period program review)Task 4F – Base period summary report (first draft)

Due date – Four months after task 3

Task Description
Modeling Methodology

Use the multi-block capability of the CPCFD code

Two situations of melt flow may require different modeling approac

The overhanging melt poses the biggest challenge

Use engineering approximation for the overhanging melt can greatproblem 

Need to know the heat and wall boundary conditions for the ovemelt and the catch basinWhat data is available from NIST for comparison? 

Block 1: Polymer object

Block 2: Overhanging melt

Block 3: Catch basin

Block 1: Polymer object

Block 2: Overhanging melt

Block 3: Catch basin

Melt drop
Surface Tension Effects

Plan to implement the continuum surface force model proposed et al (1992)

The addition of surface tension results in a pressure drop (or increfree surface and can be treated as a Dirichlet boundary condition oequation. For 2-D problem,Reference: J.U. Brackbill, D.B. Kothe, & C. Zemach, A continuum modeling surface tension, J. Comput. Phy., 100, pp335-354, 1992Need to know the surface tension coefficient for PP702 and PP652

NIST data available?If no data, how to estimate?

R p p i o

σ

− =

Pressure on the inside of the free surface

Atmosphere pressure

Radius of Surface curvatureSurface tension coefficient (temperature dependent)
 Work on Task 4A without surface tension effects Add surface tension model and complete Task 4A Complete Task 4B Complete Task 4C

Hurdles expected

Modeling overhanging meltSurface tension model

Preliminary Plan
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

Feasibility Study

(02/22/07)
Model descriptionTest case resultsSummary and future work

Outline
Study the feasibility of using CPCFD to model melt pool on catc

Case 1

catch pan not included in the computational domain assume uniform constant melt temperature at 700Kassume constant melt flow rate of 0.21 g/s (estimated from Task 3 results)assume melt flow velocity of -1.e-3 m/s perpendicular to the catch surface (estimateno heat loss to environmentcatch surface tilt angle at 0°surface tension effects neglected

Case 2

catch pan included in the computational domain but has different properties (i.e., deconductivity, and head capacity) than the polymer melt. viscosity of catch pan is set 

temperature of catch pan is 700K 

other conditions same as those in Case 1

 melt flow cumulated on the catch surface is treated as an addit

source term of the continuity equation

Additional source terms are added to momentum and energy equa

for the impacts of the cumulated melt flow

mass source term: momentum source term:energy source term: 

Model Description

m&

m(u

u) & flow

−

m(T

T) & flow

−

m& =

0.21

g /s; u flow

=

−1.e

−

3

m /s; Tflow

= 700

K
Case 1 Configuration

Initial dimension and grid

W=0.1m

H=0.0001m

L=0.05m

2-D computational domain at t = 0.1s:

L=0.05m

Grid number: 100X20

Grid lines shown here connect cell centers

Melt flow source terms are added to these two cells

free surface

no-slip wall surface

outlet BC

Y

X
Case 1 Results

t=10s

t=20s

u (m/s)

u (m/s)

Movie: 0 ~ 30 s (time interval between frames is 0.2 s)

Takes 12 hrs to run on a P5 workstation
Case 2 Configuration

Initial dimension and grid

L=0.05m

W=0.1m

2-D computational domain at t = 0.1s: Grid number: 100X30

L

Y

X

comparing to Case 1, ten additional grid direction to model the catch pan;

L X H1 uses moving grid to model cumL X H2 uses fixed grid to model catch 
Case 2 Material Properties

Viscosity (Pa-s)

Heat conductivity (W/m-K)

Density (kg/m

3)

Heat capacity (J/kg-K)

0.13 ~1.e6 1.e

0.25 1.

900 28

2400 96

Polymer melt Catch 

* data from NIST email (2/16/2007)** arbitrarily chosen to suppress the flow

Material forming catch pan has smaller heat capacity and higher heat condu

polymer melt. As a consequence, assuming 1-D heat conduction in the direcperpendicular to the catch surface inside the catch pan is invalid.

Case 2 includes the catch pan in the computational domain and may be momodeling the transient heat transfer in the catch pan.
Case 2 Results

t=10s

u (m/s)

t=20s
Case 2 Results

Movie: 0 ~ 30 s (time interval between frames is 0.2 s)

Takes 12 hrs to run on a P5 workstation
Summary and Future Work

Two cases have been successfully performed, which show that thcapable of modeling the melt flow on the catch surface under somconditionsFuture work includes:

derive heat transfer coefficient on the melt surface (free convection on horizontaincorporate more realistic boundary and initial conditions in the calculation

surface heat lossconstant temperature BC at the backside of the catch pan“cold” start (i.e., initial temperature of catch pan ~ 250°C)tilt angle greater than 0°
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

Task 4A (PP702N)

(03/21/07)
CFD modelSurface heat losses models 2-D Computational gridsTest case results

Case T1: All boundaries are adiabatic, Catch plate initial T at 600KCase T2: Free Convection on free surface and constant T (=600K) 

backside of catch plate

Case T3: Free Convection and radiative heat losses on free surfacT (=600K) at the backside of catch plateCase T4: Same as Case T3 except the tilt angle of catch plate incrto 2.5°For all cases, polymer melt mass flow rate is 0.14g/s and its tempe

Summary and future work

Outline
Constant T = 327

Catch Surface

Polymer melt (PP702N) at 0

Radiative and Convective Heat Losseson free surface

W = 7 cm

L = 9.5 in

CFD Model
Free Surface Heat Losses

Convective heat loss

Radiative heat loss

Convective heat loss:Approximate convection on the surface as “ free convection on horizontal plate” Radiative heat loss: 

q&conv

=

−

h

(Tsurf

−

Ta

), Ta

= 298

°

C

8.0

0.24

400 27

1.32

1.32

4

1 4

1

⎟ ≈⎠⎞ ⎜

⎝⎛

− ⎟ =

⎠⎞ ⎜

⎝⎛

∆

=

L

T

h

(

), 1.0, 5.67e -8

4

4

q&rad

=

−εσ

Tsurf

−

Ta

ε

=

σ

=

* Holman, J.P., Heat Transfer, McGraw-Hill, Inc., New York, 1981 (pp285)

Catch Pla

Polymer M
2-D Computational Grids

L = 0.24 m

Biased fixed grid for catch plate

Uniform (in H direction) moving grid for m

400 X 30 (L X H) grid used in Case T3 and T4

H = 0.012

L = 0.1 m

Reduced 200 X 30 (L X H) grid used in Case T1 and 
Cases Summary

Convective and radiative heat losses; Uniform T0=600K; Catch plate included; Backside constant T; NIST experiment dimension; θ = 2.5°

T4

Convective and radiative heat losses; Uniform T0=600K; Catch plate included; Backside constant T; NIST experiment dimension; θ = 0°

T3

Convective heat loss; Uniform T0=600K; Catch plate included; Backside constant T; Reduced dimension; θ = 0°

T2

03-22-2007 (cu

Adiabatic BC; Uniform T0=600K; Catch plate included; Reduced dimension; θ = 0°

T1

Adiabatic BC; Uniform T0=700K; Catch plate included; Reduced dimension; θ = 0°

2

02-23-200

Adiabatic BC; Uniform T0=700K; Catch plate not included; Reduced dimension; θ = 0°

1

Case # Description Report
Case T1

Adiabatic BC on melt free surface and catch plate surface; Uniform catch 

at 600K

Reduced computational domain (L = 0.1m) w.r.t NIST experiment (L = 0.2computational time2-dimensional grid: 200 X 30 (L X H)

0.5 inch
Case T2

Free convection heat loss on free surface included; The catch surface of plate not covered by polymer melt is adiabatic; Uniform catch plate initial Catch plate backside temperature at constant 600KReduced computational domain (L = 0.1m) w.r.t NIST experiment (L = 0.2computational time2-dimensional grid: 200 X 30 (L X H)

Convective heat loss

Adiabatic

Constant T
Case T3

Dimension of catch plate the same as NIST experimentFree convection and radiative heat loss on the free surface of the melt an

of the catch plate included; Uniform catch plate initial T at 600K; Catch pltemperature at constant 600K

2-dimensional grid: 400 X 30 (L X H)

Convective and radiative heat losses

Constant T
Case T4

Dimension of catch plate the same as experimentFree convection and radiative heat loss on the free surface and the catch 

the catch plate included; Uniform catch plate initial T at 600K; Catch plate temperature at constant 600K; Plate tilt angle 2.5°

2-dimensional grid: 400 X 30 (L X H)

Convective and radiative heat loss

Constant T

Convective and radiative heat loss

θ=2.5°
Case T4 - Snapshots 

t = 30 s (CPU = ~12 hr)

t = 300 s (CPU = ~36 hr)

t = 600 s (CPU = ~60 hr)

Calculation time steps were limited to 0.01s for t=0 ~ 30 sMaximum time steps after t = 300 s is 0.05 sSimulation was run on a 2.8 GHz Intel workstation 
Summary 

Test results from simulations using realistic boundary conditions the feasibility of the modeling methodologyMelt free surface becomes corrugated when applying heat loss There are two large peaks at the two ends of the melt in all 4 cain this presentationHeated catch surface with initial temperature above ~570K is crsuccess of the simulationsSimulation diverged when setting tilt angle to 5°Computational time steps at the beginning of the simulations (t are very small (in the order of 10-3 s). CFL = 0.1 was used in all 

and maximum time step was set to 0.05 s. It took about 12 hrs t30 s simulation and about 60 hr to complete a 600 s simulation Intel workstation. 
Future Work

Implement realistic initial temperature distribution in catch plate

Steady stateConvective and radiative heat losses at catch surfaceConstant T at backside 

Test the possibility of reducing the catch plate backside temperaturCombine simulations of the vertical resin and the catch basinSimulate PP6523 (feasibility study without surface tension model) 

Implement surface tension boundary condition on free surfaceRepeat simulations of PP702N and PP6523
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

(04/26/07)
Initial Temperature Distribution inside the catch plateCase SummaryCase T5: PP702 Simulation

Movie of polymer melt shape Temperature contour snap shots

Case T6: PP6523 Simulation

Movie of polymer melt shapeTemperature contour snap shots

Surface tension model

Surface tension force normal to the free surfaceSurface tension gradient effectsCalculation of free surface curvature

Summary and future work

Outline
Temperature Distribution of Catch Plate

Y

X

T

B

T

S

Qconv

Qrad

T

A

Qcond

1-D Steady State Energy Balance Equation:

298 K;

y 0.0127m

1.26 w/m - K; 8.0 w/m - K;

1.0; 5.67 10 w/m - K ;

( ) ( )

2 -8 2 4

4

4

=

∆

=

=

=

=

=

×

=

−

+

−

A

S

A

S

A

T

k

h

h

T

T

T

T

dy

dT

k

ε

σ

εσ

Initial catch plate temperature distribution in PP702N Simulation (T5):

( )

675 K; 588 K;

0.0127

0

6860.63

0.0127

=

=

−

≤

≤

=

− ⋅ +

T

T

y

T

T

y

B

S

B

∆

y

Initial catch plate temperature distribution in PP6523 Simulation (T6):

( )

700 K;

604 K;

0.0127

0

7570.87

0.0127

=

=

−

≤

≤

=

− ⋅ +

T

T

y

T

T

y

B

S

B
Case Summary

Convective and radiative heat losses; Uniform T0=600K; Catch plate included; Backside constant T=600K; NIST experiment dimension; θ= 2.5°; PP702N

T4

Convective and radiative heat losses; Uniform T0=600K; Catch plate included; Backside constant T=600K; NIST experiment dimension; θ= 0°; PP702N

T3

Convective heat loss; Uniform T0=600K; Backside constant T=600K; Reduced dimension; θ = 0°; PP702N

T2

03-22-200

Adiabatic BC; Uniform T0=600K; Reduced dimension; θ = 0°; PP702N

T1

Case # Description Report

Convective and radiative heat losses; Realistic initial catch plate temperature distribution; Backside constant T = 700K; NIST experiment dimension; θ = 2.5°; PP6523

T6

04-26-2007 (cu

Convective and radiative heat losses; Realistic initial catch plate temperature distribution; Backside constant T = 675K; NIST experiment dimension; θ = 2.5°; PP702N

T5
Case T5: PP702N – Shape Change

Movie time period: 2 – 360 secondsComputation time (0 – 360 secs): 15 hrs (CFL = 0.05, 

∆tmax = 0.0

Initial catch plate temperature distribution:

( )

675 K; 588 K;

0.0127

0

6860.63

0.0127

=

=

−

≤

≤

=

− ⋅ +

T

T

y

T

T

y

B

S

B
Case T5: PP702N – Temperature Snap Shots

T = 12 sec.

T5T4

T = 60 sec.

T5T4

T (K)
Case T5: PP702N – Temperature Snap Shots

T = 120 sec.

T5T4

T = 180 sec.

T5T4

T (K)
Case T5: PP702N – Temperature Snap Shots

T = 240 sec.

T5T4

T = 330 sec.

T5T4

T (K)
Case T6: PP6523 – Shape Change

Initial catch plate temperature distribution:

Movie time period: 2 – 240 secondsComputation time (0 – 240 secs) : 18 hrs (CFL = 0.05, 

∆tmax = 0.

( )

700 K;

604 K;

0.0127

0

7570.87

0.0127

=

=

−

≤

≤

=

− ⋅ +

T

T

y

T

T

y

B

S

B
Case T5: PP6523 – Temperature Snap Shots

T = 20 sec.

T6

T5

T = 60 sec.

T6

T5

T (K)
Case T5: PP6523 – Temperature Snap Shots

T = 120 sec.

T6

T5

T = 240 sec.

T6

T5

T (K)
Surface Tension as Function of Temperature - Th

2

/

3

γ

=

κ

(

TC

−

T

)

/

V

Two empirical equations for surface tension of liquids

Eotvos:

Guggenheim-Katayama 

where, 

V is the molar volume of that substance

Tc is the critical temperature 

κ is a constant for each substance.

n

T

Tc

(

1

/

)

0

γ

=

γ

−

where, 

γ0 is a constant for that substance

n is an empirical factor, whose value is 11/9 for org

http://en.wikipedia.org/wiki/Surface_tension, The critical temperature, 

material is the temperature above which distinct liquid and gas do not 
Surface Tension of Polypropylene - Experiments

0.04

T

(

C

) 29.10 (165

- 220 

C

o

o

γ

=

−

+

Schonhorn, H., and Scharpe L., “Surface tensions of Molten Polypropylen

Polymer Letters, Vol.3, 1965, p.235.
Surface Tension of Polypropylene - Experiments

0.051

T

(

C

) 30.38 (180

- 210 

C

o

o

γ

=

−

+

Kwok, D, et al., “Study on the surface tensions of polymer melts using ax

drop shape analysis”, Polymer Engineering and Science, Vol.38, 1998, p
Surface Tension of Polypropylene - Summary

The critical temperature, T

c, of a material is the temperature above 

distinct liquid and gas do not exist. 

Both Eotvos and Guggenheim-Katayama’s empirical expressions cthe fact that the surface tension reduced to zero at/above the critictemperature.. What’s is the critical temperature of PP702 and PP6523 ?

Schonhorn-Scharpe’s correlation suggests that 

T

c = 728 

o

C

Kwok’s correlation suggests that T

c = 596 

oC. 

Both Schonhorn-Scharpe and Kwok’s model are based on linear correlation, and the temperature range from 165-220 

oC 

(438 - 493 K). Eotvos’s expression suggests that linear correlation might be applicable for higher temperature (but below the critical temperature).

Phase dia

hsutesesi
Surface Tension CFD Mode – Surface Tension Fo

Surface Tension Forces (in the normal direction of the free surface

Young-Laplace:Implementation in CFD as a body forceSurface tension 

γ as a function of temperature (Thomas J. Ohlemiller 1/16/07 em

Model only applies to cells forming the free surface

⎟⎟⎠⎞

⎜⎜⎝⎛

∆

=

+

1 2

1 1

R

R

P

γ

∆P is the pressure difference over the inteR1, R

2 are the principal radii of curvature 

γ is the surface tension 

nn

n

n

Fvol

P

A

n

A

n

vv

v

=

∇

•

=

=

−

∆ ⋅ ⋅ =

− ⋅ ⋅ ⋅

ˆ; ˆ

ˆ

ˆ

κ

γ

κ

κ is the curvature of the fre

is the normal vector of th

A is the surface area

n

v

0 040 2910 (165 - 220 C)

0.051 30.38 (180 - 210 C)

=

−

+

°

=

−

+

°

. T .

T

γγ Philips HNZ020 polypropyl

Eastman Epolene D-10 atac
Surface Tension CFD Mode – Surface Tension Gr

Surface Tension Gradient Effects (in the tangential direction of the 

Shear stress at the surface resulting from the variation of surface tension*:Implementation in CFD as a momentum source Model only applies to cells forming the free surface

τ γ

τ

γ

μ

=

≡

∂

∂

dd

n

U μ is the viscosity; 

U is the tangential velocity componentn is the normal direction of the free suτ is the tangential direction of the free 

τ

τ

τ

γ

γ

γ

τ

γ

γ

τ

τ

τ

τ

τ

vv

v

v

v

=

= ⋅

= ⋅

= ⋅ =

+

ˆ

ˆ

S

V

S

V

S i j

v y

u x

x y

S

u , Sv are the additional source terms in the u- and v-momentum equation

is the tangential vector of the free surface

V is the cell volume

τ

v

* W.A. Sirignano and I. Glassman, Combustion Sci. and Tech., Vol. 1, pp307-312, 
Surface Tension CFD Mode – Calculate Surface C

Define the free surface passing through (x1,y1), (x

2,y

2), (x

3,y

3) using parabo

y=g(x)=a1+a

2x+a

3x2. 

Solve linear algebraic equations for coefficients

ai

y1=a1+a

2

x1+a

3x1

2.

y2=a1+a

2

x2+a

3x2

2.

y3=a1+a

2

x3+a

3x3

2.

Define a functional f(x,y)=y-g(x)=y- a1- a

2x- a

3

x

2, the normal to the surface 

The curvature of the free surface may be expressed as:

2 1/ 2 (1 ( ) )

( )

| | ˆ

g

x

j

g

x i

f

f

+ ′

− ′ =

∇

∇

=

r

r

n

2 3/ 2 (1 ) ˆ

x

xx

g

g

+

=

∇⋅ =

n

κ

(x1,y1) (x2,y2) (x3,y3)

n

y=f(x)=a1+a2x+a3x2

x

y

Reference:http://web.mit.edu/1.63/www/le
Summary 

Realistic initial condition of catch plate has been implemented in modelStarting from the realistic initial condition, two test cases (T5 anbeen performed for PP702N and PP6523. Comparing the resulta previous case (T4, see 3/22/07 report) that started from a uniftemperature initial condition, the following observations can be m

 In order for the CFD code to run with realistic initial condition, the b

temperature must be increased from 600 K (used in T4) to 675 K aT5 and T6, respectively

 Comparing to case T4, the melt flow in T5 and T6 expands much f

the higher backside temperatures, and

 The free surfaces are much smoother in T5 and T6 than that in T4

 Because of the different viscosity-temperature relationships, the pr

behaviors of the two polymer melts are different in terms of free su

and expansion speeds

 Surface tension model has been developed
Future Work

Implement surface tension model into CPCFDRepeat test cases T5 and T6 with surface tension effectsCombined simulations of the vertical resin and the catch basin for PP702N and PP6523
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

(05/30/07)
Implement Surface Tension ModelCase SummaryStudy the Impact of Time Step Size

Snap shots of free surface shape at certain time pointsMean expansion velocities vs. time 

Study the Impact of Surface Tension Effects

Movie of polymer melt shapeTemperature contour snap shotsMean expansion velocities vs. timeMomentum source terms due to surface tension effects

Summary and future work

Outline
Surface Tension CFD Mode – Surface Tension Fo

Surface Tension Forces (in the normal direction of the free surface

Young-Laplace:Implementation in CFD as a body forceSurface tension 

γ as a function of temperature (Thomas J. Ohlemiller 1/16/07 em

Model only applies to cells forming the free surface

⎟⎟⎠⎞

⎜⎜⎝⎛

∆

=

+

1 2

1 1

R

R

P

γ

∆P is the pressure difference over the inteR1, R

2 are the principal radii of curvature 

γ is the surface tension 

nn

n

n

Fvol

P

A

n

A

n

vv

v

=

∇

•

=

=

−

∆ ⋅ ⋅ =

− ⋅ ⋅ ⋅

ˆ; ˆ

ˆ

ˆ

κ

γ

κ

κ is the curvature of the fre

is the normal vector of th

A is the surface area

n

v

γ

=

−

0.040

T

+ 29.10 (165 - 220

°C) Eastman Epolene D-10 at
Surface Tension CFD Mode – Surface Tension Gr

Surface Tension Gradient Effects (in the tangential direction of the 

Shear stress at the surface resulting from the variation of surface tension*:Implementation in CFD as a momentum source Model only applies to cells forming the free surface

τ γ

τ

γ

μ

=

≡

∂

∂

dd

n

U μ is the viscosity; 

U is the tangential velocity componentn is the normal direction of the free suτ is the tangential direction of the free 

τ

τ

τ

γ

γ

γ

τ

γ

γ

τ

τ

τ

τ

τ

vv

v

v

v

=

= ⋅

= ⋅

= ⋅ =

+

ˆ

ˆ

S

A

S

A

S i j

v y

u x

x y

S

u , Sv are the additional source terms in the u- and v-momentum equation

is the tangential vector of the free surface

V is the cell volume

τ

v

* W.A. Sirignano and I. Glassman, Combustion Sci. and Tech., Vol. 1, pp307-312, 
Calculate n1 the normal vector to line passing through (x1,y1), (x2,y2), and c

normal vector to line passing through (x2,y2), (x3,y3). n can be computed fro

Then the curvature at node point (x2, y2) can be approximated as. Similarly, we can obtain the curvature at node point (x3,y3), (x4,y4).... If we 

curvature at a cell center, we may just simply take the average of curvature points. 

y

n n

x

n n

y

n

x

n

y

n

x

n

x x y y

x y x y

∆

−

+

∆

−

=

∆

∆

+

∆

∆

=

∂

∂

+

∂

∂

= ∇⋅ =

2 1 2 1

2 n κ

Surface Tension CFD Model – Calculate Surface C

κ2

(x1,y1)

(x2,y2)

(x3,y3)

x y

n1 n2

(xc2,yc2) (xc3,yc3)

κ3

κ4

(xc4,yc4)(x4,y4)

(x,y) : node position(xc,yc) : cell position

∆x

∆y
Case Summary

∆tmax = 0.02 Included Inclu

CFL = 0.5

T12 Same as T11

Not Included Not Inc

∆tmax = 0.02

CFL = 0.5

Convective and radiative heat losses; Backside constant T = 700K; θ = 2.5°; 

T11

∆tmax = 0.02 Included Not Inc

CFL = 0.5

T10 Same as T7

∆tmax = 0.02 Included Inclu

CFL = 0.5

T9 Same as T7

Not incNot inc

Not includedNot included

∆tmax = 0.1 

CFL = 0.5

Convective and radiative heat losses; Backside constant T = 675K; θ = 2.5°; 

T7

∆tmax = 0.02 

CFL = 0.5

T8 Same as T7

Surface Tension 

effect

(normal direction)

Surface T

effe

(tangential 

Cases Boundary and Initial Conditions Time Step 

T = T − 6860.63⋅( ) y + 0.0127 B

T = T − 7570.87⋅(y + 0.0127) B

16.030 - 250T11

CPU time (hr) 8.5 16.5 21.0 21.0 Time Period (second) 30 - 250 30 - 250 30 - 250 30 - 250 Cases T7 T8 T9 T10 
Case T7 vs T8: Impact of Time Step

T7

T8

T = 100 sec.

T8

T7

T = 200 sec.

T7

T = 300 sec.

CPU time (30 – 300 sec):

Case T7 – 8.5 hrCase T8 – 19.5 hr
Case T7 vs T8: Impact of Time Step

0.0E+00

1.0E-042.0E-043.0E-044.0E-045.0E-04

50.00 100.00 150.00 200.00 250.00 3

Time (s)

Mean Velocity of Melt Front (m/s)T7, dtmax = 0.1, Right Front T7, dtmax = 0.1, Left Front T8, dtmax = 0.02, Right Front T8, dtmax = 0.02, Left Front
Case T8 vs T9: Impact of Surface Tension Effects

Case T8, 30 – 250 sec, CPU time = 16.5 hr

Case T9, 30 – 250 sec, CPU time = 21 hr

PP
Case T8 vs T9: Impact of Surface Tension Effects

0.0E+00

2.0E-044.0E-046.0E-048.0E-04

50.00 100.00 150.00 200.00 25

Time (s)

Mean Velocity of Melt Front (m/s)

T8, Right FrontT8, Left FrontT9, Right FrontT9, Left Front

PP
Case T8 vs T9: Snap Shots of Temperature Conto

T8

T9

t = 150 sec

t = 250 sec

T8

T9
Case T9: Source Terms due to Surface Tension F

n

n

Fvol

P

A

v

=

∇

•

=

=

−

∆ ⋅ 

ˆ; ˆ

κ

Fv

t = 50 sec
Case T9: Source Terms due to Surface Tension G

τ

τ

τ

γ

γ

γ

τ

τ

τ

vv

v

=

= 

= 

= 

ˆ

SSS 

v 

u 

S

T (K)
Case T9 vs. T10: Effects of Surface Tension Grad

T1

T9

t = 250 sec

The magnitude of the surface tension gradient source terms are msmaller than that of the surface tension force source terms The impact of neglecting the surface tension gradient effects is neunder the simulated condition.
Case T11 vs T12: Impact of Surface Tension Effe

Case T11, 30 – 250 sec, CPU time = 16.0 hr

Case T12, 30 – 250 sec, CPU time = 40 hr 

PP
Case T11 vs T12: Impact of Surface Tension Effe

PP

0.0E+00

2.0E-044.0E-046.0E-048.0E-04

50.00 100.00 150.00 200.00 250.00

Time (s)

Mean Velocity of Melt Front (m/s)

T11, Right FrontT11, Left FrontT12, Right FrontT12, Left Front
Cases T11 & T12: PP6523 – Temperature Snap Sh

t = 160 sec. T11

T12

t = 250 sec. 

T11T12

T12

T11
Case T12: PP6523N – Surface Tension Source Te

Force in the Normal Direction t = 160 sec.

Shear Force in the Tangential Direction
Case T9 vs T12: Expansion Rates Comparison

0.0E+00

2.0E-044.0E-046.0E-048.0E-04

50.00 100.00 150.00 200.00 2

Time (s)

Mean Velocity of Melt Front (m/s)

T9, PP702N, Right FrontT9, PP702N, Left FrontT12, PP6523, Right FrontT12, PP6523, Left Front
Summary 

Surface tension model has been implemented into CPCFDA parametric study on time step size has been performed

reducing maximum time step size from 0.1s to 0.02s increases the mean ex

by 2% at 50 s and by 10% at 300 s.

reducing maximum time step size from 0.1s to 0.02s doubles the run-time

The inclusion of surface tension effects

makes the free surface much smootherincreases the expansion velocities (moving downward) by 30% for PP702N 

PP6523 at the end of the calculations (i.e., t = 250 sec)

The melt front (moving downward) of PP6523 moves faster than under the simulated conditionsThe surface tension gradient effect owning to surface temperatuis negligible comparing to the surface tension force in the normathe free surface under the simulated condition.
Viscosity vs. Temperature

1.00E-01

1.00E+001.00E+011.00E+021.00E+031.00E+041.00E+051.00E+06

50 150 250 350 450 550 650

Temperature (K)

Viscosity (Pa-s)

PP702NPP6523

1.00E-01

1.00E+001.00E+011.00E+021.00E+031.00E+04

575 600 

Tempera

Viscosity (Pa-s)

Free surface temperature varies in the range of 575 K to 650 KViscosity of PP6523 is much larger (by one order of magnitudeof PP702N from 575 K to about 645 KFor temperature greater than 645 K, viscosity of PP6523 becosmaller than that of PP702N
Future Work

Parametric study of grid resolutionCombined simulations of the vertical resin and the catch basin for PP702N and PP6523Simulation of gasification process without flow?Finalize source code and documentationFinal report
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

(06/27/07)
Time Step Sensitivity Study

Snap shots of free surface shape at certain time pointsMean expansion velocities vs. time

Grid Sensitivity Study

Snap shots of free surface shape at certain time pointsMean expansion velocities vs. time 

Summary and future work

Outline
Case Summary

T19 ∆tmax = 0.02; CFL = 0.5 533X30 15.0 Not IncluT18 ∆tmax = 0.045; CFL = 0.5 400X30 15.0 Not inclu

48.036.021.021.017.034.015.0

8.5

T9 ∆tmax = 0.02; CFL = 0.5 400X30 Include

T16 ∆tmax = 0.02; CFL = 0.5 600X30 IncludeT17 ∆tmax = 0.02; CFL = 0.5 800X30 IncludeT15 ∆tmax = 0.02; CFL = 0.5 800X30 Not IncluT14 ∆tmax = 0.02; CFL = 0.5 600X30 Not IncluT13 ∆tmax = 0.01; CFL = 0.5 400X30 Not Inclu

Not incluNot inclu

400X30

T7 ∆tmax = 0.1; CFL = 0.5 400X30T8 ∆tmax = 0.02; CFL = 0.5

Cases Time Step Grid Surface Tensi

CPU time (hr)*

* estimated, physical time period from 30 s to 250 s

** all calculations are for PP702N
Impact of Time Step on Free Surface Shape

T7 T8

T = 50 sec. T8 T7

T = 125 sec.

T7

T = 250 sec.

T13

T13
Impact of Time Step on Spreading Velocity

0.0E+00

1.0E-042.0E-043.0E-044.0E-045.0E-04

50.00 100.00 150.00 200.00 

Time (s)

Mean Velocity of Melt Front (m/s)

T7, dtmax = 0.1, Right FrontT7, dtmax = 0.1, Left FrontT8, dtmax = 0.02, Right FrontT8, dtmax = 0.02, Left FrontT13, dtmax = 0.01, Right FrontT13, dtmax =0.01, Left Front
2.80E-042.85E-042.90E-042.95E-043.00E-043.05E-043.10E-043.15E-043.20E-04

0 0.002 0.004 0.006 0.008 0.01 

Squre of Maximum Time Step (∆tmax

2

)

Mean Velocity of Right Front at 250 sec. (m/s)

Impact of Time Step on Spreading Velocity
Impact of Grid Size on Free Surface Shape

T8

T15

T = 100 sec.

T = 200 sec.

T14

Surface Tension Effects not included
0.0E+00

5.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04

50 100 150 20

Time (s)

Mean Velocity (m/s) 400X30, Right Front 400X30, Left Front

600X30, Right Front600X30, Left Front800X30, Right Front800X30, Left Front

Impact of Grid Size on Spreading Velocity

Surface Tension Effects not included
Impact of Grid Size on Free Surface Shape

T9

T17

T = 50 sec.

T17

T = 130 sec.

Surface Tension Effects included

T16
Impact of Grid Size on Spreading Velocity

Surface Tension Effects included

0.0E+00

5.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-044.0E-04

50 60 70 80 90 100 110 120 

Time (s)

Mean Velocity of Melt Front (m/s) 400X30, Right Front 400X30, Left Front 600X30, Right Front 600X30, Left Front 800X30, Right Front 800X30, Left Front
Impact of Grid Size on Spreading Velocity

1.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04

0 0.2 0.4 0.6 0.8 1 1.2

Square of Normalized Grid Size (∆X

2

)

Mean Right Front Velocity at 130 sec. (m/s) Surface Tension Effects not Included Surface Tension Effects Included

“Accurate solution” by Extrapolation
Summary 

A parametric study on time step size has been performed

reducing maximum time step size increases the mean expansion velocity reducing maximum time step size from 0.1s to 0.02s doubles the run-time

A grid sensitivity study has been performed

Fluctuations of free surface become stronger as grid size decreasesThe inclusion of surface tension effects suppresses the fluctuationsThe inclusion of surface tension effects leads to smaller differences among fvelocities using different grids
Future Work

Combined simulations of the vertical resin and the catch basin for PP702N and PP6523Simulation of gasification process without flow?Finalize source code and documentationFinal report
A Computational Model For 

Fire Growth and Spread On 

Thermoplastic Objects

REACTION

INTERNATIONAL

ENGINEERING

Task 4: Include Melt Pool on Catch Surf

(07/20/07)
Task 4A – model PP702 

Time dependent melt flow rate obtained from Task3 (Case 3R, see 12/05/06 rCatch plate tilt angle = 2.5°

 Task 4B – model PP6523

 Time dependent melt flow rate obtained from Task3 (Case 3I, see 09/26/06 re

 Catch plate tilt angle = 2.5°

 Task 4C – model PP702 with a different tilt angle

 Melt flow rate the same as Task 4A

 Catch plate tilt angle = 0.0°

 Summary and future work

Outline
Case Summary

* Surface tension model included in all simulations

4C PP702N From Case 3R 0.01 s 400X30 0.0°

4B PP6523 From Case 3I 0.02 s 400X30 2.5°

4A PP702N From Case 3R 0.02 s 400X30 2.5°

Case Polymer Melt flow rate Time step Grid Tile angle
0

0.05

0.1

0.15

0.2

0.25

0.3

0 100 200 300 400 5

Time (s

Melt flow rate (g/s)

Task 4A – Melt flow rate

0

50

100150200250300

0 200 400 600 800

Time (s)

Mass Loss (g)

Mass loss history 

- solid line = total mass loss- dashed line = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasification

Melt flow rate calculated as the first derivative of the dashed line

Case 3R
Task 4A – Melt flow on catch plate

Physical time period: 200 sec – 460 sec CPU time: 19 hrs

• The small melt flow rate onto the catch plate at the beginning ( < 100 sec. ) caused ca

• melt temperature became too low due to heat loss and small flow rate

• To overcome the difficulty, the melt flow rate on to the catch plot from 0 sec. to 200 se

modeled using a constant melt flow rate to create an initial amount of mass on the ca

1) Compute the melt flow using a large constant melt flow rate (i.e., 0.14 g/s) for 30 seconds. Thi

about 5 grams of melt on the catch plate.

2) Stop simulation 3) Restart simulation using the time dependent melt flow rate profile, starting at 200 sec. 

4) Continue simulation until melt front approaches end of plate
Task 4A - Spreading velocity

0.0E+00

1.0E-042.0E-043.0E-044.0E-045.0E-04

200.00 250.00 300.00 350.00 400.00 

Time (s)

Mean Velocity of Melt Front (m/s)

4A, Right Front4A, Left Front
0

50

100150200250300

0 500 1000 1500 2000

Time Exposed (s)

Mass Loss (g)

Task 4B – Melt flow rate

Mass loss history 

- solid line = total mass loss- dashed line = mass loss occurring as melt flow- difference between solid and dashed line = mass loss due to gasification

Melt flow rate calculated as the first derivative of the dashed lineNote: The melt flow rate in Case 3I is one order of magnitude smaller than that in Case 3R

Case 3I

0

0.0020.0040.0060.008

0.01

0.0120.0140.0160.018

0 400 800 

Time 

Melt Flow Rate (g/s)
Task 4B – Melt flow on catch plate

Physical time period: 700 sec – 1400 sec CPU time: 36 hrs

• The small melt flow rate at the beginning ( < 500 sec. ) caused case to crash

• melt temperature became too low due to heat loss and small flow rate

To overcome the difficulty, the melt flow rate on to the catch plot from 0 sec. to 700 sec

modeled using a constant melt flow rate to create an initial amount of mass on the ca

1) Compute the melt flow using a large constant melt flow rate (i.e., 0.14 g/s) for 30 seconds. Thi

about 5 grams of melt on the catch plate.

2) Stop simulation 3) Restart simulation using the time dependent melt flow rate profile, starting at 700 sec. 

4) Continue simulation until melt front approaches end of plate
Task 4B - Spreading velocity

0.0E+00

2.0E-054.0E-056.0E-058.0E-051.0E-041.2E-041.4E-04

700 800 900 1000 1100 1200 1300 

Time (s)

Mean Velocity of Moving Fronts (m/s) Right Front Left Front
Task 4C – Melt flow on catch plate

Physical time period: 200 sec – 460 sec CPU time: 40 hrs

• Melt flow rate the same as that in Task 4A• Used same start-up procedure as Task 4A. 
Task 4C - Spreading velocity

0.0E+00

5.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04

200 300 400 500 

Time (s)

Mean Velocity of Moving Fronts (m/s) Right Front Left Front

Note: identical curves for right front and left front
Summary and Future Work

Required Simulations in Task4 have been performed

Code encountered problem when started with time dependent melt flow rate is too small at the early stage, which leads to low melt temperature due to heA restart methodology was developed to overcome the problem. We first run using a relatively large constant melt rate for a certain time period. We then 

simulation with the time dependent melt rate. Mass conservation is used to d

restart point of the time dependent melt rate

Future Work

Finalize source code and documentationFinal report
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 1 

APPENDIX B 

CPCFD Code Documentation 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 2 

1.0 Instructions for using the CPCFD code to solve problems of Task 1, 2 and 3

1) Initial Grid Generation 

 

 Before running CPCFD, the user must edit, compile and run mesh.f to generate the initial 

computational grid. The grid is stored in two input files GRID and INLET. Currently, the code 

is written for the Task 1 problem and the grid is 19X200 (ghost cells not included) for the demo 

case. 

2) Input File 

 File input.dat stores all the parameters that control the simulation. The format of a typical 

input.dat for Task 1 is listed below: 

 T T F T F F F !sol_u v w p tk ep om 

 0 0 0 !nrs nmf itcm 

 1 1 1 1 1 1 !NSWPu p tk ep om scalars 

 0.85 1.00 0.8 0.8 0.3 1.00 !RELAXu p tk ep den scalars 

 3 1 1 1 3 !ISCMu tk ep om scalars 

 1.013e5 0 20 1 0 !p_ref i,j,k,ib_init 

 0.0 -9.8 0.0 !gx, gy, gz 

 2000 2.5 !itermax sormax 

 0 1000 1 1 !restart nitout nitflow nitsca 

 1 !is_unsteady 

 400.0 2.d-0 1.0d0 1 !t_end dt t_outp time_level 

 1.d-1 !CFL number 

 2.4d3 2.98d2 9.0d2 2.5d-1 !Cp, temp_init, density, heat conductivity 

 30000.d0 1.25d6 !incident heat flux, gasification heat 

 2 !ipolymer, 1 - PP23K; 2 - PP702N; 3 - PP6523 

 The first line contains logical flags that determine which equations are to be solved. For 2-D 

polymer melting problem, the user must not change this line. For 3-D problem, the third 

parameter should be changed to "T" that enables the code to solve the "w velocity". 

 The second line contains parameters related to reactive flow simulation, and should be kept the 

same when modeling polymer melting problems. 

 

 The third line contains the number of TDMA sweeps used in solving the linear equations. 

 The fourth line contains the under-relaxation factors used in solving different equations. The 

first under-relaxation factor corresponds to the momentum equation, which may need 

adjustments during calculation to avoid divergence or to speed up convergence. This can be 

done by changing the value in the ACTIVE file. The other relevant under-relaxation factors are 

the second one and the last one, corresponding to the pressure equations and the temperature 

equation, respectively. Currently they are set to be 1.0 (i.e., no under-relaxation), and have not 

caused any problem in all the test runs. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 3 

 The fifth line controls the differencing schemes for the equations. 1 – upwind scheme; 2 -- 

central differencing scheme; 3 -- mixing scheme; 4 -- exponent scheme. 

 The sixth line contains the reference pressure (Pa), and the index of the reference grid point. 

 The seventh line contains the components of the gravity vector 

 The eighth line contains the maximum number of outer iteration in each time step (itermax) 

and the residual limit (sormax) used to determine convergence. For example, if sormax=1.0, 

then the solutions must converge to a level that at least six significant digits remain unchanged 

from iteration to iteration. Whereas if sormax = 6, only one significant digit remain unchanged. 

 The ninth line contains, from left to right, the flag for restart, the iteration interval for 

diagnostic output, the number of inner iterations for solving the momentum and pressure 

equations, and the number of inner iterations for solving the temperature equation. If the restart 

flag is 1, then the code will try to read in the restart files at the beginning of the simulation. If it 

is 0, the simulation will start from scratch. For restart, see section 8) for detail. 

 The tenth line contains the flag for simulating unsteady problem, 1 for time accurate 

unsteady problem and 0 for steady problem. For Task 1, the flag should remain 1. 

 The eleventh line contains, from left to right, the ending time of the simulation, the maximum 

time step allowed, the time interval for solution output, and the time marching scheme (0 for 

first order backward Euler scheme, and 1 for second order 3-level scheme). 

 The twelfth line contains the CFL number used to determine time step. Larger CFL number 

leads to longer time step. 

 The last two lines assigns case parameters: specific heat, initial temperature, density, heat 

conductivity, and incident heat flux, gasification heat and a flag for choosing polymer. 

3) ACTIVE File 

 Once the simulation starts, an ACTIVE file will be written out, and can be used to control the 

simulation in the fly. A typical ACTIVE file reads 

 

 1000.0 !End time 

 2000 !Max iteration # 

 1 !Max flow iteration # 

 1 !Max temperature iteration # 

 0.3 !density under-relaxation factor 

 100 !number of outer iterations between diagnostic screen output 

 0.85 1.0 !under-relaxation factors for momentum and pressure equations 

 1.0 !under-relaxation factor for temperature equation 

 2.0 !maximum time step 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 4 

 0.2 !CFL number 

 1.0 !output time interval 

 2.5 !residual upper limit 

 

 See comments at the end of each line for explanation. 

 

 The minimum number of global iterations (default value) in each time step is 10. 

 

4) Case Dependent Source Files 

 The Fortran source code includes the following files: 

 

 cpcfd.f -- main program 

 fmb3.f -- containing all driving subroutines 

 fblock.f -- implementations of models and solver. 

 A Makefile can be used to compile and link the source code. The user may need to modify the 

Makefile based on the available operating system and Fortran compiler. At REI, the code is 

compiled under Linux system using a Portland Group Fortran Compiler (pgf90). 

 Certain parts of the source code are case dependent, and need modification when setting up 

different cases. They are contained in fblock.f and are listed below. For more detail, see the code. 

 

 subroutine fblck_setup --- set material properties such as k, Cp, density, and initial 

 temperature. 

 subroutine fblck_mvg --- implement re-meshing scheme 

 subroutine fblck_pcorrct --- implement cell face mass flux correction due to moving 

 mesh and boundary 

 subroutine fblck_bc_sca --- apply heat flux boundary condition 

 function vislam --- implement viscosity-temperature relation 

 In the demo, the code is configured to mode PP702 exposed to 30,000 W/m2 constant heat flux 

with radiative and convective heat loss. 

5) Run the Simulation from Scratch 

 Following are the steps to run the demo simulation from scratch 

 a) Create two sub-directories (./FV and ./output) under working directory. 

 b) Compile and run mesh.f. Copy GRID and INLET to the working directory. 

 c) Compile the source code and copy the executable CPCFD.exe to working directory. 

 d) Edit and Copy input.dat to working directory. Make sure the restart flag is turned off. 

 e) launch CPCFD.exe 

 Typical screen output looks like the following: 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 5 

 ................ 

 ................ 

 ................ 

 Time = 0.20868E+03 W(y=0) = 0.28049E-01 W(y=0.25) = 0.99674E-02 

 Time = 0.20868E+03 totv = 0.58634E-03 Mass Loss Ratio = 0.61879E-01 Total mass loss 

rate = 0.21649E+00 Melt flow fraction = 0.82360E+00 

 mass conservation: 0.34031181E-03 0.35697476E-08 0.76451802E-09 0.13137932E-03 

0.28340087E-04 0.50003122E-03 

 0.51208481E-02 0.29436213E-09 0.51208481E-01 0.29436213E-09 

 TIME STEP: 4800 time = 208.69372 dt = 0.00512 

 TIME STEP: 4800 time = 208.69372 dt = 0.00512 

 ITER = 1 Iter_flow = 1 RESIDUAL: M, U, V, W 

 1.09237 0.17468 0.69568 -22.96195 

 TIME STEP: 4800 time = 208.69372 dt = 0.00512 

 ITER = 1 Iter_scalar = 1 Residual: SCALARS 

 1.62199 

 MUVS 1 1.09237 0.17468 0.69568 1.62199 

 MUVS 10 0.96373 0.29299 0.70663 -4.84325 

 Time = 0.20869E+03 W(y=0) = 0.28049E-01 W(y=0.25) = 0.99669E-02 

 Time = 0.20869E+03 totv = 0.58634E-03 Mass Loss Ratio = 0.61893E-01 Total mass loss 

rate = 0.21649E+00 Melt flow fraction = 0.82360E+00 

 mass conservation: 0.34030748E-03 0.35694999E-08 0.76450220E-09 0.13138289E-03 

0.28340851E-04 0.50003122E-03 

 0.51201375E-02 0.29437275E-09 0.51201375E-01 0.29437275E-09 

 ................ 

 ................ 

 ................ 

 The first line shows, from left to right, the current time, object thickness at y=0 (bottom) and at 

y=0.25 (top). 

 The second line shows, from left to right, the current time, volume of the calculation domain, 

mass loss ratio, total mass loss rate and melt flow fraction during the current time step. 

 The third line shows information regarding mass conservation. From left to right, they are 

volume remaining in the domain, volume flowing out during the previous time step, volume 

gasifying into the environment, cumulated volume flowing out of the domain, cumulated volume 

gasifying into the environment, and the sum of the remaining volume and the cumulated volumes 

leaving the domain. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 6 

 The fourth line shows information regarding time step selection. From left to right, they are 

actual time step used, maximum time step based on diffusion, maximum time step based on 

convection, and maximum time step based on both diffusion and convection. 

 The fifth line shows the current time step, current time, and time step size. 

 The sixth to eighth lines show the residuals of the continuity and momentum equations (listed 

in the eighth line) after one global iteration. 

 The ninth to eleventh lines show the temperature equation residual. 

 The interval of global iterations for residual output can be adjusted in the ACTIVE file 

 The next two lines are summary of the residuals for each solved equations. 

6) Output Data 

 

 The solutions that can be visualized by graphic tools are written to the ./FV directory. The 

data are in PLOT3D format. The data are written out every n seconds, where n is the output 

time interval defined in input.dat file and can be adjusted in the ACTIVE file. The files are 

organized as follows: 

 f3dmbXXX.fun --- stores the solution 

 f3dmbXXX.grd --- stores the multi-block grid 

 f3dmb.nam --- stores the name of the solved quantities 

 In the file name, "XXX" represents the time in seconds. For example, f3dmb500.fun contains 

the solution at t=500 sec. In addition to these data files, file f3dmb0.fun and f3dmb0.grd always 

store the latest solution after each time step, which can be used to monitor the solution without 

having to wait the solution to reach an specific time point. 

 

 When the simulation completes, the final solution is written to the following files: f3dmb.fun 

and f3dmb.grd 

 In REI, FieldView is used to visualize the data. When reading the data, choose FORMATTED 

for FILE FORMAT; choose 3-D for COORDINATES; and choose MULTI-GRID for DATA 

FORMAT. 

7) Other Diagnostic Output 

 Several files are written to the working directory to monitor the simulation process. 

 File surf.dat contains information about the free surface location. The data are organized into 

three columns. The first column stores the elevation of the grid points on the free surface, and the 

second and third columns store the corresponding coordinates in the horizontal direction at the 

beginning and at the end of the current time step. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 7 

 File massloss.dat contains information about mass conservation, which are also shown on the 

screen. See section 5) for detail. 

8) Restart 

 Restart file are stored in the ./output directory. After each output time point, two files are 

written out, they are 

 GRIDXXX and RESTRTXXX 

 where "XXX" represents the time in seconds. 

 To restart, the user must copy these files to the working directory and rename them as GRID 

and RESTRT, respectively. For example, if the user wants to restart the simulation at t=140 sec. 

the files GRID140 and RESTRT140 must be copied to the working directory and renamed to 

GRID and RESTRT. 

 Before restart, make sure the restart flag in input.dat is set to 1. 

 Delete GRID and RESTRT if the simulation is intended to start from scratch. In this case, 

GRID must be generated by executing the initial mesh generator. See section 1) for detail. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 8 

2.0 Instructions for using the CPCFD code to solve problems of Task 4

1) Initial Grid Generation 

 

 Before running CPCFD, the user must edit, compile and run mesh.f to generate the initial 

computational grid. The grid is stored in two input files GRID and INLET. For the Task 4 demo 

case, there are 30 grid points in the X direction and 400 grid points in the Y direction (ghost cells 

not included). The catch plate consists of 10X400 fixed grid points. The polymer melt has 

20X400 grid points and these points are allowed to move in the X direction. The initial thickness 

of the polymer melt is set to be 1.e-4 m. For more details about the computational grid, see final 

reports. 

2) Input File 

 File input.dat stores all the parameters that control the Task 4 simulation. The format of a 

typical input.dat for Task 4 is listed below: 

 T T F T F F F !sol_u v w p tk ep om 

 0 0 0 !nrs nmf itcm 

 1 1 1 1 1 1 !NSWPu p tk ep om scalars 

 0.85 1.00 0.8 0.8 0.3 1.00 !RELAXu p tk ep den scalars 

 3 1 1 1 3 !ISCMu tk ep om scalars 

 1.013e5 0 20 1 0 !p_ref i,j,k,ib_init 

 0.0 -9.8 0.0 !gx, gy, gz 

 2000 2.5 !itermax sormax 

 0 1000 1 1 !restart nitout nitflow nitsca 

 1 !is_unsteady 

 400.0 2.d-0 1.0d0 1 !t_end dt t_outp time_level 

 1.d-1 !CFL number 

 2.4d3 2.98d2 9.0d2 2.5d-1 !Cp, temp_init, density, heat conductivity 

 30000.d0 1.25d6 !incident heat flux, gasification heat 

 2 !ipolymer, 1 - PP23K; 2 - PP702N; 3 - PP6523 

 -1.0 650.0 !mass flow rate, melt temperature 

 1 !surface tension model 

 The first line contains logical flags that determine which equations are to be solved. For 2-D 

polymer melting problem, the user must not change this line. For 3-D problem, the third 

parameter should be changed to "T" that enables the code to solve the "w velocity". 

 The second line contains parameters related to reactive flow simulation, and should be kept the 

same when modeling polymer melting problems. 

 

 The third line contains the number of TDMA sweeps used in solving the linear equations. 

 The fourth line contains the under-relaxation factors used in solving different equations. The 

first under-relaxation factor corresponds to the momentum equation, which may need 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 9 

adjustments during calculation to avoid divergence or to speed up convergence. This can be 

done by changing the value in the ACTIVE file. The other relevant under-relaxation factors are 

the second one and the last one, corresponding to the pressure equations and the temperature 

equation, respectively. Currently they are set to be 1.0 (i.e., no under-relaxation), and have not 

caused any problem in all the test runs. 

 The fifth line controls the differencing schemes for the equations. 1 -- upwind scheme; 2 -- 

central differencing scheme; 3 -- mixing scheme; 4 -- exponent scheme. 

 The sixth line contains the reference pressure (Pa), and the index of the reference grid point. 

 The seventh line contains the components of the gravity vector - gx, gy, and gz. The catch 

plate tilt angle beta is also defined by the gravity vector as tan(beta) = gy/gx. 

 The eighth line contains the maximum number of outer iteration in each time step (itermax) 

and the residual limit (sormax) used to determine convergence. For example, if sormax=1.0, 

then the solutions must converge to a level that at least six significant digits remain unchanged 

from iteration to iteration. Whereas if sormax = 6, only one significant digit remain unchanged. 

 The ninth line contains, from left to right, the flag for restart, the iteration interval for 

diagnostic output, the number of inner iterations for solving the momentum and pressure 

equations, and the number of inner iterations for solving the temperature equation. If the restart 

flag is 1, then the code will try to read in the restart files at the beginning of the simulation. If it 

is 0, the simulation will start from scratch. For restart, see section 8) for detail. 

 The tenth line contains the flag for simulating unsteady problem, 1 for time accurate 

unsteady problem and 0 for steady problem. For Task 1, the flag should remain 1. 

 The eleventh line contains, from left to right, the ending time of the simulation, the maximum 

time step allowed, the time interval for solution output, and the time marching scheme (0 for 

first order backward Euler scheme, and 1 for second order 3-level scheme). 

 The twelfth line contains the CFL number used to determine time step. Larger CFL number 

leads to longer time step. 

 The thirteenth to fifteenth lines assign case parameters: specific heat, initial temperature, 

density, heat conductivity, and incident heat flux, gasification heat and a flag for choosing 

polymer. It should be noted that, when simulating melt flow on the catch plate, the initial 

temperature in the data file is only used to set the backside temperature of the catch plate. Initial 

temperature distribution within the catch plate is hardwired in the code currently (see lines 794 - 

800 of fblock.f) and the final report 

 The sixteenth line assigns a constant mass flow rate of the melt (kg/s) and the melt 

temperature (K). To use time dependent mass flow rate in the calculation, the user must do the 

following: 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 10 

 A. set the mass flow rate in the input file to -1.0 

 B. create an ASCII data file "massrate.dat" in the following format 

 

 11 -------- Total number of table entry 

 0.0 7.81E-05 -------- Time (s) and mass flow rate (kg/s) 

 5.01332 8.16E-05 

 10.0187 8.60E-05 

 15.02898 9.20E-05 

 20.03494 9.84E-05 

 25.03807 1.04E-04 

 30.04405 1.07E-04 

 35.04669 1.08E-04 

 40.05655 1.08E-04 

 45.06237 1.09E-04 

 50.07059 1.10E-04 

 The code will tabulate the table contained in "massrate.dat" to determine the time dependent 

mass flow rate 

 The last line is a flag for the surface tension model. To apply the surface tension model, set the 

flag to 1. Otherwise, set it to 0. 

 

3) ACTIVE File 

 Once the simulation starts, an ACTIVE file will be written out, and can be used to control the 

simulation in the fly. A typical ACTIVE file reads 

 1000.0 !End time 

 2000 !Max iteration # 

 1 !Max flow iteration # 

 1 !Max temperature iteration # 

 0.3 !density under-relaxation factor 

 100 !number of outer iterations between diagnostic screen output 

 0.85 1.0 !under-relaxation factors for momentum and pressure equations 

 1.0 !under-relaxation factor for temperature equation 

 2.0 !maximum time step 

 0.2 !CFL number 

 1.0 !output time interval 

 2.5 !residual upper limit 

 

 See comments at the end of each line for explanation. 

 

 The minimum number of global iterations (default value) in each time step is 10. 

 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 11 

4) Case Dependent Source Files 

 The Fortran source code includes the following files: 

 

 cpcfd.f -- main program 

 fmb3.f -- containing all driving subroutines 

 fblock.f -- implementations of models and solver. 

 A Makefile can be used to compile and link the source code. The user may need to modify the 

Makefile based on the available operating system and Fortran compiler. At REI, the code is 

compiled under Linux system using a Portland Group Fortran Compiler (pgf90). 

 Certain parts of the source code are case dependent, and need modification when setting up 

different cases. They are contained in fblock.f and are listed below. For more detail, see the code. 

 

 subroutine fblck_setup --- set material properties such as k, Cp, density, and initial 

 temperature. 

 subroutine fblck_mvg --- implement re-meshing scheme 

 subroutine fblck_pcorrct --- implement cell face mass flux correction due to moving 

 mesh and boundary 

 subroutine fblck_bc_sca --- apply heat flux boundary condition 

 function vislam --- implement viscosity-temperature relation 

 subroutine SurfCurveA --- calculates free surface curvature 

 The surface tension model is implemented in fblck_bc_sca, fblck_calcu, and fblck_calcv 

 

 In the demo, the code is configured to mode the Task 4 problem and can not be used for Task 1 

to 3. 

 

5) Run the Simulation from Scratch 

 Following are the steps to run the demo simulation from scratch 

 a) Create two sub-directories (./FV and ./output) under working directory. 

 b) Compile and run mesh.f. Copy GRID and INLET to the working directory. 

 c) Compile the source code and copy the executable CPCFD.exe to working directory. 

 d) Edit and Copy input.dat to working directory. Set t_end to 30.0, mass flow rate to 1.4e-4. 

Make sure the restart flag and the surface tension model flag are turned off. 

 e) Launch CPCFD.exe 

 Once the simulation of the first 30 seconds completed, do the following 

 

 a) Copy GRID29 and RESTRT29 from ./output to the working directory and rename them as 

GRID and RESTRT 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 12 

 b) Open input.dat. Set t_end to 400.0 and mass flow rate to -1.0. Turn on the restart flag and 

the surface tension model flag. Save input.dat. 

 c) Copy massrate.dat to the working directory. 

 c) Launch CPCFD.exe again. 

 Typical screen output looks like the following: 

 ................ 

 ................ 

 ................ 

 Time = 0.32455E+02 W(j=2) = 0.10000E-03 W(j=200) = 0.45568E-02 VR = 0.14790E-03 

VL = 0.16639E-03 

 melt flow rate = 1.0748186883438162E-004 

 0.20000000E-01 0.35000003E-09 0.99929105E-01 0.35000003E-09 

 

 TIME STEP: 149 time = 32.49467 dt = 0.02000 

 TIME STEP: 149 time = 32.49467 dt = 0.02000 

 ITER = 1 Iter_flow = 1 RESIDUAL: M, U, V, W 

 2.42088 1.73416 2.11876 -22.40645 

 TIME STEP: 149 time = 32.49467 dt = 0.02000 

 ITER = 1 Iter_scalar = 1 Residual: SCALARS 

 1.43313 

 MUVS 1 2.42088 1.73416 2.11876 1.43313 

 MUVS 10 2.49579 1.77480 2.13870 0.31786 

 Time = 0.32475E+02 W(j=2) = 0.10000E-03 W(j=200) = 0.45544E-02 VR = 0.14781E-03 

VL = 0.16628E-03 

 melt flow rate = 1.0748586672349618E-004 

 0.20000000E-01 0.35000003E-09 0.10109939E+00 0.35000003E-09 

 ................ 

 ................ 

 ................ 

 The first line shows, from left to right, the current time, the thickness of the melt layer at the 

edge of the catch plate, the thickness of the melt layer at the center of the catch plate, the moving 

velocity of the right front and the moving velocity of the left front. 

 The second line shows the melt flow rate at the current time step. 

 The third line shows information regarding time step selection. From left to right, they are 

actual time step used, maximum time step based on diffusion, maximum time step based on 

convection, and maximum time step based on both diffusion and convection. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 13 

 The fourth line shows the current time step, current time, and time step size. 

 The fifth to seventh lines show the residuals of the continuity and momentum equations 

(listed in the seventh line) after one global iteration. 

 The eighth to tenth lines show the temperature equation residual. 

 The interval of global iterations for residual output can be adjusted in the ACTIVE file 

 The next two lines are summary of the residuals for each solved equations. 

6) Output Data 

 

 The solutions that can be visualized by graphic tools are written to the ./FV directory. The data 

are in PLOT3D format. The data are written out every n seconds, where n is the output time 

interval defined in input.dat file and can be adjusted in the ACTIVE file. The files are 

organized as follows: 

 f3dmbXXX.fun --- stores the solution 

 f3dmbXXX.grd --- stores the multi-block grid 

 f3dmb.nam --- stores the name of the solved quantities 

 In the file name, "XXX" represents the time in seconds. For example, f3dmb500.fun contains 

the solution at t=500 sec. In addition to these data files, file f3dmb0.fun and f3dmb0.grd always 

store the latest solution after each time step, which can be used to monitor the solution without 

having to wait the solution to reach an specific time point. 

 

 When the simulation is completed, the final solution is written to the following files: 

f3dmb.fun and f3dmb.grd 

 In REI, FieldView is used to visualize the data. When reading the data, choose FORMATTED 

for FILE FORMAT; choose 3-D for COORDINATES; and choose MULTI-GRID for DATA 

FORMAT. 

7) Other Diagnostic Output 

 File surf.dat contains information about the free surface location. The data are organized into 

three columns. The first column stores the elevation of the grid points on the free surface, and the 

second and third columns store the corresponding coordinates in the horizontal direction at the 

beginning and at the end of the current time step. 

8) Restart 

 Restart file are stored in the ./output directory. After each output time point, two files are 

written out, they are 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 14 

 GRIDXXX and RESTRTXXX 

 where "XXX" represents the time in seconds. 

 To restart, the user must copy these files to the working directory and rename them as GRID 

and RESTRT, respectively. For example, if the user wants to restart the simulation at t=140 sec. 

the files GRID140 and RESTRT140 must be copied to the working directory and renamed to 

GRID and RESTRT. 

 Before restart, make sure the restart flag in input.dat is set to 1. 

 Delete GRID and RESTRT if the simulation is intended to start from scratch. In this case, 

GRID must be generated by executing the initial mesh generator. See section 1) for detail. 
REI 4714 A Computational Model For Fire Growth & Spread On Thermoplastic Objects B - 15 

3.0 Instructions for adjusting control parameters to avoid divergence

It is possible that, during the course of a simulation, the code diverges. To solve the problem, it is 

always helpful to make sure that all the restart files are properly stored. If a divergence happens, 

one can always restart the case from a previous time step after adjusting some of the control 

parameters in the input data file. The following steps can be taken to make the code less prone to 

divergence: 

1. reducing the maximum time step size 

2. reducing the CFL number that also controls the time step size 

3. reducing the under-relaxation factors used in solving the momentum, pressure, and 

energy equations 

If the problem persists, try to restart from an even earlier time step and repeat the above steps. 

The proper values of these control parameters are case dependent. It is necessary to adjust these 

numbers on a trial-and-error base. 

The control parameters can also be adjusted in the ACTIVE file without terminating the case. By 

monitoring the screen output or the output log file, it is often possible to spot the on-set of 

divergence. For example, if the number of iterations to reach convergence in a time step 

increases significantly compared to the number of iterations required in previous time steps, it 

may indicate that the simulation has become unstable and divergence may occur within the next 

couple time steps. To avoid the divergence, try the three steps mentioned above in the ACTIVE 

file. 

Sometimes, divergence is caused by factors such as using extreme boundary conditions and thus 

divergence can not be solved by simply adjusting the control parameters. For example, in the 

simulations performed for Task 4 if the backside temperature of the catch plate is below a certain 

value, the case will diverge due to the worsening grid quality (i.e., grid skewness) developed 

during the simulation. If this occurs, one has to study the exact cause of the problem and fix it.