materials

Article
In Situ, Real-Time Temperature Mapping and Thermal FE
Simulations of Large-Format 3D Printed PETG/CF Vertical Wall
Felipe Robles Poblete 1 , Matthew Ireland 1 , Lucinda Slattery 2 , William G. Davids 3
and Roberto A. Lopez-Anido 3, *

1 Advanced Structures and Composites Center (ASCC), University of Maine, Orono, ME 04469, USA;
[Link]@[Link] (F.R.P.); [Link]@[Link] (M.I.)
2 Department of Physics and Astronomy, University of Maine, Orono, ME 04469, USA;
[Link]@[Link]
3 Department of Civil and Environmental Engineering, University of Maine, Orono, ME 04469, USA;
[Link]@[Link]
* Correspondence: rla@[Link]; Tel.: +1-207-581-2119

Abstract: This work focuses on simulating the thermal history of a vertical wall consisting of a
thermoplastic composite material, poly(ethylene terephthalate) glycol (PETG) with short carbon
fiber reinforcement, manufactured using a Big Area Additive Manufacturing (BAAM) system. The
incremental deposition process used in additive manufacturing, which corresponds to the repeated
deposition of hot material onto cooler material, contributes to the presence of residual stresses and
part warping. The prediction of these mechanisms is dependent on thermal history of the part, and
the major motivation of this work was to improve the accuracy of finite element (FE) models used to
quantify the thermal history of large-format additively manufactured parts. Thermocouples were
placed throughout the part at varying heights to measure temperature as a function of time. The
FE model developed found a thermal contact conductance between the printed part and the bed of
10 W/m2 K and convection coefficient values that linearly varied from 3 to 15 W/m2 K through the
wall height when making a temperature comparison with the output from the thermocouples. It is
also demonstrated that the FE model with a constant convection coefficient under-predicts model
Citation: Robles Poblete, F.; Ireland,
temperature at the beginning of the manufacturing process when compared against the model with
M.; Slattery, L.; Davids, W.G.;
a variable convection coefficient. The impact of this difference was seen in the stress values, which
Lopez-Anido, R.A. In Situ, Real-Time
were larger for the model with a constant convection coefficient. Finally, a correlation equation was
Temperature Mapping and Thermal
FE Simulations of Large-Format 3D
derived which allows the findings to be generalized to other vertical structures manufactured on
Printed PETG/CF Vertical Wall. the BAAM. In summary, this work offers valuable insights on material characterization, real-time
Materials 2023, 16, 6486. https:// thermocouple placement, and FE modeling of large-format additively manufactured parts.
[Link]/10.3390/ma16196486
Keywords: additive manufacturing; finite element; modeling; conductance; convection; additive
Academic Editors: Xiaochao Liu
manufacturing; thermoplastic; polymer
and Lei Shi

Received: 15 August 2023

Revised: 31 August 2023
Accepted: 4 September 2023 1. Introduction
Published: 29 September 2023
Large-format extrusion-based additive manufacturing is a technology that, in recent
years, has become widespread in the fabrication of composite applications in the marine
and construction industries, allowing for greater design flexibility while reducing lead times
Copyright: © 2023 by the authors.
and costs [1–3]. The Big Area Additive Manufacturing (BAAM) [4] system developed at
Licensee MDPI, Basel, Switzerland. Oak Ridge National Laboratory in collaboration with Cincinnati Inc. has been successfully
This article is an open access article used to manufacture large parts with a variety of thermoplastic composites [5,6]. Both
distributed under the terms and the BAAM system and desktop-scale Fused Filament Fabrication (FFF) extrude heated
conditions of the Creative Commons thermoplastic material along programmed tool paths to manufacture parts on a layer-by-
Attribution (CC BY) license (https:// layer basis [6]. Unlike FFF, BAAM uses a single-screw extruder to melt polymer pellets
[Link]/licenses/by/ and force the molten material through a nozzle via a pressure differential [7] instead of
4.0/). resistively melting a thin filament feedstock. Single-screw extruders enable the use of

Materials 2023, 16, 6486. [Link] [Link]

Materials 2023, 16, 6486 2 of 18

thermoplastic materials at a relatively lower cost and at faster deposition rates with mass
throughputs up to 50 kg/h. With regard to the architecture of manufactured parts, BAAM-
and FFF-produced components are similar, although the former produces parts which are
an order of magnitude larger with bead dimensions at or above 10 mm. An example of a
large-format additively manufactured structure is BioHome3D, which is a 56 m2 modular
house manufactured using a recyclable biopolymer filled with wood fiber [8].
Performance objectives for 3D-printed parts create demand for materials to exhibit
functionalities [9] including improved electrical and thermal conductivity, mechanical
strength, and stiffness at relatively low cost [5,10]. To accommodate this demand, re-
searchers have attempted to mix different types of fillers, such as metal [10], glass fibers [11],
and vapor grown short carbon fibers [11,12] into the polymer matrix. Although the
macrostructure of large 3D-printed parts alone contributes to thermal and mechanical
anisotropy due to layer-wise deposition [12,13], this phenomenon is accentuated in short-
fiber composites. Fibers with varying aspect ratios tend to align in the print direction,
significantly affecting homogenized material properties [7]. Moreover, fiber alignment has
also been seen to vary within the printed bead itself [14].
For a fiber-reinforced thermoplastic polymer, the cooling behavior of the deposited
material is governed by heat transfer to the environment due to convective and radiative
heat losses as well as conduction between beads and layers [15]. The rate of cooling governs
both the phase change from viscoelastic fluid to solid and inter-bead bond quality [16]. The
combination of these processes impacts the formation of residual stresses and deformations
within the part [17], affecting the shape of the extrudate [18] and subsequent mechanical
properties [19,20]. High-quality characterization of the associated thermal history is there-
fore required to ensure robust prediction of outcomes from the manufacturing process.
Moreover, new additive manufacturing technologies such as 4D printing rely on the pro-
gramming of different extrudate temperatures throughout printing to govern the polymer
structure shape transformation over time [21,22].
Experimentally obtained thermal history data are often captured and reported by
means of infrared (IR) radiation thermography [14,23–25]. The initialization of an IR
camera for data capture, which is a necessary component of thermographic measurement
for accurate temperature reporting, requires information about the scene and subject to
correctly correlate as-measured radiance with as-reported temperatures. This requisite
information includes an emissivity parameter, which varies according to the material
surface roughness, the temperature dependence of the material’s emissive response, the
angle of incidence between the subject surface normal and camera optical axis, and the
line-of-sight distance between subject and camera. Similarly, knowledge of the scene
temperature is required during initialization to accurately account for the proportion of as-
measured radiance due to reflection from the subject [26]. These factors are often neglected
or simplified, as in the case of a constant emissivity value, the practice of which imposes
error on the temperature data reported by IR cameras [25]. By contrast, thermocouples
are commonly used in research and industry to measure temperature with relatively
simple sources of inaccuracy and across different processes. Previously, thermocouples
have been embedded in small-scale additively manufactured parts for in situ temperature
characterization [27,28].
Predictive tools that incorporate the coupled impacts of bed temperature, ambient
temperature, and material properties on thermal history are also necessary. Layer-by-
layer deposition models have been developed ranging from simple axisymmetric 1D
transient heat transfer models [14,29] to 3D finite element models (FE) [30,31]. The finite
difference method has also been used to numerically model temperature variation for
the FFF processes, including large-format additive manufacturing, due to reduced com-
putation costs when compared against FE implementations [32–34]. Recently, a coupled
thermo-mechanical numerical model to determine a suitable combination of the parame-
ters that avoids the collapse of the deposited layer under self-weight was developed [35].
FE-based methods have represented the thermal history of 3D-printed parts based on
Materials 2023, 16, 6486 3 of 18

element activation [36], accounting for temperature-dependent material properties [37]

and presented features that allow modeling heat transfer at time scales small enough
to capture rapid cooling events [38]. In particular, the commercial FE software Abaqus
([Link] accessed on 1 Au-
gust 2023) with additive manufacturing capabilities has been used to model complex
3D-printed parts, such as cellular structures with homogenized material properties [39]
and thin-walled tubular structures [40]. Abaqus has also been applied to model the ther-
mal history, final deformed shape, and residual stresses in additively manufactured parts
comprised of acrylonitrile butadiene styrene (ABS) polymer [41], ABS with short carbon
fibers [15,42], polyphenylene sulfide (PPS) polymer with carbon fibers [43,44], and met-
als [45]. The majority of published research on the topic of thermal modeling for FFF rely
on the use of a constant convection coefficient.
This work combines in situ temperature measurements obtained from an additively
manufactured part with candidate FE models of the manufacturing process. Candidate
models were compared against experimental data, and the FE implementation that mini-
mized error was found to require a non-constant convection coefficient in order to accurately
capture the thermal history of the part. Finite element analysis was used to model the
complete thermal history of a large-format 3D printed vertical wall made of poly(ethylene
terephthalate) glycol (PETG) with short carbon fiber (CF) reinforcement. PETG is recog-
nized for its manufacturability with glass transition and melting temperatures of 85 ◦ C and
260 ◦ C, respectively [46,47], qualifying the material as a good candidate for thermal and
structural characterization. The accuracy of the thermal model was enhanced by real-time
temperature data gathered by thermocouples embedded in the part during the manufactur-
ing process. The temperature correlation between experimentally obtained and numerically
generated data facilitated the characterization of conductance between the part and print
bed, as well as convective heat transfer between the part and the environment, comprising
process model features which were found to substantially impact the development of resid-
ual stresses. Finally, a correlation equation was derived based on the analysis of the wall
manufactured with PETG/CF material and tested on a separate wall manufactured with
ABS/CF. The necessity of this study is driven by the tendency for large-scale additively
manufactured parts to fracture and/or develop significant distortion during manufacturing
due to the accumulation of residual stresses [48–51]. Hence, the objective of this work is to
improve the accuracy of FE models intended to capture the thermal behavior of large-scale
polymer AM during fabrication via in situ temperature measurements.

2. Materials and Methods

2.1. Printing Process Information
Part manufacturing was executed on the BAAM machine stationed in the Advanced
Structures and Composites Center at the University of Maine campus in Orono, Maine.
A prismatic vertical wall was chosen for geometric simplicity and to facilitate parametric
convection studies via the measurement of temperature variations along the height. The
wall was manufactured with Techmer Electrafil 1711 PETG, which is compounded with
18% carbon fiber volume fraction. The magnitude of fiber volume fraction was not chosen
to satisfy any specific criteria, but it is typical for materials provided by the supplier, and
prior studies have used products with similar amounts [14,25]. The average carbon fiber
length and diameter were 163 µm and 7 µm, respectively. The wall consisted of a single
bead with the first layer extended laterally to form a brim for improved stability. The
initial manufacturing process parameters were based on the layer time utilized in a prior
publication with similar geometry [14] and modified to mitigate the overall deformation and
debonding between layers. The wall dimensions and manufacturing process parameters
are given in Table 1.
Type K thermocouples were manually installed between layers of the wall to capture
the temperature at the interfaces. Thermal history was obtained below the first layer at the
part/bed interface and at layers 38, 77, 116, 155, and 194. A plywood scaffolding structure
Materials 2023, 16, 6486 4 of 18

was designed to hold the thermocouple leads and utilized to prevent forces due to gravity
from pulling thermocouples out of position during solidification of the extrudate.

Table 1. PETG/CF wall dimensions and printing process information.

Dimensions Value Printing Process Information Value

Wall Height (m) 1 Deposition Temperature (◦ C) 200
Wall Length (m) 0.75 Bed Set Temperature (◦ C) 90
Bead Width (mm) 15.875 Dryer Temperature (◦ C) 60
Bead Height (mm) 5.08 Feed Rate (m/min) 3.24
Bead Count 1 Deposition/Layer Time (s) <168/13.9
Layer Count 197 Manufacturing Time (hours) <10

2.2. Interlayer Thermocouples—Final Position Measurement

Although the approximate locations of embedded thermocouples were known from
visual inspection, accurate positional measurements were taken to verify the quality of
contact with the extrudate. A Quantum Max FaroArm® (Faro Technologies, Lake Mary,
FL, USA) was used to provide a detailed 3D scan of each wall for X–Z (length, height)
location determination. A coordinate system was chosen to denote the locations of the
thermocouples, and the planar geometry of the walls was leveraged accordingly. The
coordinate system origin in X, Y and Z was chosen to be the first identifiable point where
extrusion begins, the midpoint of the wall thickness, and the center point at the part/bed
interface, respectively. Thermocouple locations are reported in reference to this coordinate
system and were used for comparison with model data.
End-mill removal of the as-manufactured material revealed bond quality with the
surrounding polymer as well as the relative location of the thermocouple within the bead.
Overall, 30% of the thermocouples embedded in the PETG/CF wall exhibited some aspect
of poor-quality bonds (C3, AI0, and AI3) as determined by visual inspection of the contact
between the thermocouple lead wires and the extrudate. The set of thermocouples observed
to have good bonds with their surrounding polymer was used as sources for comparison
with model data.
Figure 1 shows the position coordinates together with their calculated uncertainty
values and labels for each interlayer thermocouple in the PETG/CF wall with respect to
the direction that material was deposited in a given layer. Location uncertainty was char-
acterized by disruptions in the external surface of the extrudate, which created regions of
scan data devoid of information due to occlusion. Positions were determined by averaging
the extreme values of the disruption in X, Y, and Z orientations. Finally, the average was
then subtracted from the maximum value to determine the associated uncertainty. Thermo-
couples labeled as C0 and C1 were placed at the part/bed interface and are not included in
Figure 1. The letter difference in thermocouple labels (C5 ofand
Materials 2023, 16, x FOR PEER REVIEW 19
AI) denotes sampling rates of
1 and 2 Hz, respectively.

Figure 1. Measured positions of each interlayer thermocouples on the PETG/CF wall.

Figure 1. Measured positions of each interlayer thermocouples on the PETG/CF wall.
2.3. PETG/CF Material Characterization
The material characterization procedures utilized in this work were based on a pro‐
posed roadmap for testing the same type of additively manufactured short‐fiber compo‐
site materials [52]. Thermomechanical and mechanical property data were obtained as in‐
puts to the FE models. Material property data were generated from test specimens exer‐
cised from a different part manufactured in PETG/CF with the same deposition tempera‐
ture profile, deposition speed, and nominal bead dimensions as the wall print. Aligning
Materials 2023, 16, 6486 5 of 18

2.3. PETG/CF Material Characterization

The material characterization procedures utilized in this work were based on a pro-
posed roadmap for testing the same type of additively manufactured short-fiber composite
materials [52]. Thermomechanical and mechanical property data were obtained as inputs to
the FE models. Material property data were generated from test specimens exercised from
a different part manufactured in PETG/CF with the same deposition temperature profile,
deposition speed, and nominal bead dimensions as the wall print. Aligning the processing
conditions in parts manufactured for material property characterization with the processing
conditions employed for experimental prints is intended to control for uncharacterized
process effects.
Density measurements were performed with the Specific Gravity Method according
to the ASTM D792-20 [53]. In total, there were 54 samples (18 samples cut in each X, Y,
and Z orientation). Each sample was 17 × 17 × 7 mm3 with the short axis parallel to the
orientation of interest. The final density parameter utilized in the FE model was the sample
set average value of ρ = 1271.185 kg/m3 .
Specific heat (Cp ) measurements were performed with a TA Differential Scanning
Calorimeter instrument (DSC2500—TA Instruments, New Castle, DE, USA) according to
ASTM D3418-21 [54]. Five samples weighing at least 5 mg were tested to verify consistency
in the measured response. Table 2 shows a subset of the values utilized in the FE model,
which was an average of the five samples from 25 to 225 ◦ C.

Table 2. Average of specific heat out of five different PETG/CF samples for x, y, and z orientations.

Temperature (◦ C) 25 60 75 90 125 200 225

Cp (J/kg◦ C) 763.8 903.2 1154.4 1181.6 1268.8 1442.4 1496.6

Thermal conductivity at room temperature was determined using the transient plane
source (TPS) method according to ISO 22007-2 [55] which utilizes direct thermal diffusivity
measurement. In total, 81 paired combinations out of 54 samples (18 samples for each
X, Y, and Z orientations) were tested at room temperature of approximately 25 ◦ C. Each
sample was cut with the dimensions of 17 × 17 × 7 mm3 . Capturing the temperature
dependence of the conductivity response for the same range of temperatures adopted
for specific heat characterization is ideal, however, limitations to equipment functionality
prevented this level of fidelity. After averaging the measurement data among all samples for
each orientation, the orthotropic thermal conductivity values at room temperature adopted
in the FE model were 0.59, 0.48. and 0.35 W/m2 K for the X, Y, and Z directions, respectively.
Coefficient of thermal expansion (CTE) values were obtained by using a TA Ther-
momechanical Analyzer (TMA Q400—TA Instruments, New Castle, DE, USA) according
to ASTM E831-19 [56]. In total, 5 samples for each orientation, X (4.9 × 4.9 × 8.2 mm3 ),
Y (4.9 × 4.9 × 6.9 mm3 ), and Z (4.9 × 4.9 × 4.1 mm3 ), were tested. Strain measurements
from the test were preserved for temperatures below the glass transition temperature
(Tg = 74.4 ◦ C) of the material, which was determined according to ASTM D7028-07 [57].
At temperatures above Tg , thermally-induced strains were assumed to be constant. The
CTE for each orientation was obtained by dividing the strain measurement values by the
difference between a temperature of interest and the reference temperature of Tref = 20 ◦ C.
The average strain curves and their derived CTE curves for each orientation are shown
in Figure 2.
Elastic response in X and Z orientations as a function of temperature was measured by
a TA Dynamic Mechanical Analyzer (DMA850—TA Instruments, New Castle, DE, USA)
according to ASTM D5023-15 [58]. Three rectangular specimens (49 × 2 × 10.6 mm3 ) for X
and Z orientations were tested in flexure as a beam. The elastic response in Y orientation
was assumed to be the same as the Z orientation response for simplification. Values for shear
moduli and Poisson’s ratios were obtained from published tensile and compressive test
data [59,60]. The subsequent room temperature orthotropic elastic response was used as a
 Coefficient of thermal expansion (CTE) values were obtained by using a TA Thermo‐
mechanical Analyzer (TMA Q400—TA Instruments, New Castle, DE, USA) according to
ASTM E831‐19 [56]. In total, 5 samples for each orientation, X (4.9 × 4.9 × 8.2 mm3), Y (4.9
× 4.9 × 6.9 mm3), and Z (4.9 × 4.9 × 4.1 mm3), were tested. Strain measurements from the
test were preserved for temperatures below the glass transition temperature (Tg = 74.4 °C)
Materials 2023, 16, 6486 6 of 18
of the material, which was determined according to ASTM D7028‐07 [57]. At temperatures
above Tg, thermally‐induced strains were assumed to be constant. The CTE for each ori‐
entation was obtained by dividing the strain measurement values by the difference be‐
reference definition for the multi-factor approach [61] in order to represent the temperature
tween a temperature of interest and the reference temperature of Tref = 20 °C. The average
dependence of the elastic stiffness, as shown in Table 3. The elastic stiffness values were
strain curves and their derived CTE curves for each orientation are shown in Figure 2.
assumed to be constant for temperatures at and above 74.2 ◦ C.

Figure 2. Average
Figure 2. Averageof
oftotal
totalCTE
CTEand
andstrain
strainvalues
valuesout
outofof5 5different
differentPETG/CF
PETG/CFsamples
samplesfor
forx,x,y,y,and
and z
orientations.
z orientations.

TableElastic response
3. Orthotropic in properties
elastic X and Z orientations astemperature
as a function of a function used
of temperature was measured
in the FE model.
by a TA Dynamic Mechanical Analyzer (DMA850—TA Instruments, New Castle, DE,
Temperature (◦ C) Exx (MPa) USA)Eyy
according
(MPa) toEzz
ASTM
(MPa) D5023‐15 [58]. Three
Gxy (MPa) rectangular
Gxz (MPa) specimens
Gyz (MPa) νxy(49 ×ν2 yz × 10.6
νxzmm3)
20 12,100 ± 15for X and
2720Z orientations were tested
2720 ± 10 1110 in flexure1110
as a beam. The
1134 elastic0.38
response
0.39 in Y orien‐
0.38
32.8 6876 ± 18 tation was
2612 assumed to be the same as the Z orientation response for simplification. Values
2612 ± 12 2509 2509 972 0.32 0.34 0.32
44.3 6761 ± 19 for shear
2569
moduli and ± 12
2569 Poisson’s 2467 were obtained
ratios 2467 955
from published 0.32
tensile0.34 0.32
and compres‐
54.2 6640 ± 24 2523 2523 ± 13 2423 2423 938 0.32 0.34 0.32
64.5 6389 ± 25 sive test data [59,60].
2427 2427The
± 13 subsequent
2332 room temperature
2332 orthotropic
903 elastic0.34
0.32 response
0.32 was
74.2 4289 ± 25 used as1630
a reference definition
1630 ± 13 for1565
the multi‐factor
1565 approach606[61] in order
0.32 to0.34represent
0.32 the
temperature dependence of the elastic stiffness, as shown in Table 3. The elastic stiffness
values wereFE
2.4. Thermal assumed
Model oftothe
bePETG/CF
constantWall
for temperatures at and above 74.2 °C.
Thermal models of the single-bead PETG/CF wall manufactured on the BAAM were
Table 3. Orthotropic elastic properties as a function of temperature used in the FE model.
implemented in Abaqus/CAE 2021.HF8. The models utilized the Additive Manufactur-
Temperature (°C) Exx (MPa) Eyy (MPa)
ing (AM) moduleEzzof(MPa) Gxy (MPa)
Abaqus that Gxz (MPa)
drives sequential Gyz (MPa)
element activationνby
xy means νyz of an
νxz
20 event series.
12,100 ± 15 2720 An in-house
2720 ± 10MATLAB code was used
1110 1110to generate
1134the event series0.39
0.38 from the
0.38
32.8 G-Code-based
6876 ± 18 2612 definition
2612 ± 12of the toolpath
2509 given to the BAAM numerical
2509 972 controller.
0.32 0.34While
0.32
the cross-section of the extruded layers is approximately elliptical, for simplicity, the models
44.3 6761 ± 19 2569 2569 ± 12 2467 2467 955 0.32 0.34 0.32
assume a rectangular bead cross-section.
54.2 6640 ± 24 2523 2523 ± 13 2423 2423 938 0.32 0.34 0.32
The wall was meshed with linear hexahedral heat transfer elements (DC3D8) with
64.5 6389 ± 25 2427 2427 ± 13 2332
a seed interval set equal to the layer height of 5.076 mm. The same interval 0.34
2332 903 0.32 was used0.32
74.2 4289 ± 25 1630 length
for element 1630
and± 13 1565
width, producing 1565comprised606
a mesh of cubic0.32
elements.0.34 Mesh
0.32
convergence studies were performed separately to ensure the mesh density chosen for
this analysis was acceptable. The bed was modeled using the same DC3D8 heat transfer
element in direct contact with the brim. The thickness for the bed geometry was 1.6 mm
corresponding to the thickness of the ABS sheet placed on the bed for printing. A density of
1140 kg/m3 , a thermal conductivity of 0.17 W/m2 K, and a specific heat of 1640 J/kg-K [62]
were used for the ABS/CF sheet. Figure 3 shows the FE model including the wall and bed
along with an image of the wall mesh.
 gence studies were performed separately to ensure the mesh density chosen for this ana
ysis was acceptable. The bed was modeled using the same DC3D8 heat transfer elemen
in direct contact with the brim. The thickness for the bed geometry was 1.6 mm corre
sponding to the thickness of the ABS sheet placed on the bed for printing. A density o
1140 kg/m3, a thermal conductivity of 0.17 W/m2K, and a specific heat of 1640 J/kg‐K [62
Materials 2023, 16, 6486 7 of 18
were used for the ABS/CF sheet. Figure 3 shows the FE model including the wall and be
along with an image of the wall mesh.

Figure 3. Image Figure

of the FE modelofofthe
3. Image theFE
PETG/CF
model of wall on top ofwall
the PETG/CF the on
BAAM
top ofbed
the (left)
BAAM and
bedan(left)
imageand an imag
of the mesh usedofon
the mesh
the used onwall
PETG/CF the PETG/CF
(right). wall (right).

Thermal analysis Thermal

usedanalysis used a heat
a heat transfer steptransfer
with a step
fixed with
timea fixed time increment
increment of 10 seconds
of 10 s. This
time-step valueThiswastime‐step value wasand
tested separately tested separately
selected and it
because selected
provides because it provides
a balance betweena balance be
computation time tween
andcomputation time and
solution accuracy. Thesolution accuracy.
top surface of theThe
bedtop wassurface of thea bed
assigned was assigne
convec-
tion coefficientaofconvection 2 K, whichofwas
2.55 W/mcoefficient W/m2K, which
2.55estimated for a was estimated
horizontal for asurface
planar horizontal planar sur
[63],
and an emissivityface value
[63], and an emissivity
of 0.92 [64]. Ambientvaluetemperature
of 0.92 [64]. Ambient temperature
for the convective for the convectiv
coefficient
applied to the coefficient
ABS/CF applied
bed wastocaptured
the ABS/CF bed thermocouple
from was captured from data thermocouple data at the steady
at the steady-state
regime. A fixed temperature boundary condition of 74.5 C, as measured by thermocouplesby thermo
state regime. A fixed temperature boundary ◦ condition of 74.5 °C, as measured
installed at the couples
part/bed installed at the
interface, waspart/bed interface,
set for the bottom was
andsetside
for the bottom
surfaces of and side surfaces
the bed, and of th
bed, and the same temperature was used as an
the same temperature was used as an initial condition for the entire bed. initial condition for the entire bed.
The thermal history data from FE models were extracted at selected nodes whic
The thermal history data from FE models were extracted at selected nodes which
correspond to the measured locations of the thermocouples. Synchronization in time wa
correspond to the measured locations of the thermocouples. Synchronization in time was
necessary to accurately compare the experimental and model‐generated values. Data ob
necessary to accurately compare the experimental and model-generated values. Data
tained from thermocouples exhibit a “ramp onset” feature, which is defined as the mo
obtained from thermocouples exhibit a “ramp onset” feature, which is defined as the
ment in time when the extrudate is deposited over the thermocouple and a sharp rise i
moment in time when the extrudate is deposited over the thermocouple and a sharp rise
temperature is observed. This feature of the thermocouple time series was aligned wit
in temperatureitsisequivalent
[Link] Thisactivation
feature offeature
the thermocouple
in the FE model. time series
After was aligned the inter
synchronization,
with its equivalent
polation of FE data was conducted such that the number of sample points the
nodal activation feature in the FE model. After synchronization, was equal i
interpolation ofpreparation
FE data was forconducted such that
root‐mean‐square the analysis.
(RMS) number of sample points
Conductance was equal
was initially varied in th
in preparation for root-mean-square (RMS) analysis. Conductance was initially varied in
the FE models, and the temperature results were compared to the thermocouple (TC) data
obtained at the part/bed interface (C0 and C1). The convection study was then carried
out by comparing FE model temperature data with temperature data for all subsequent
thermocouples that exhibited good quality contact with the extrudate (C2, C4, C5, C6, C7,
AI1, and AI2).
In the following section, results are presented on how conductance and convection
values in the FE model were found from fitting the model with experimental data. Each
fitting was assessed based on RMS analysis; final conductance and convection values
minimized the error between experimentally obtained thermal histories and their model-
predicted equivalents. Several RMS time window sizes from 5 to 300 s were tested, issuing
similar results. A visualization of the study progression is shown in Figure 4.
 and AI2).
In the following section, results are presented on how conductance and convection
values in the FE model were found from fitting the model with experimental data. Each
fitting was assessed based on RMS analysis; final conductance and convection values min
imized the error between experimentally obtained thermal histories and their model‐pre
Materials 2023, 16, 6486 8 of 18
dicted equivalents. Several RMS time window sizes from 5 to 300 seconds were tested
issuing similar results. A visualization of the study progression is shown in Figure 4.

Figure
Figure 4. Diagram 4. Diagram
of the of the plan
plan followed followed
during during conductance
conductance andstudy
and convection convection
of thestudy of the FE models
FE models.

3. Results [Link]
Results and Discussion
3.1. Wall/Bed Thermal Conductance
3.1. Wall/Bed ThermalStudy
Conductance Study
For the conductance
For the study,
conductancethe convection
study, thecoefficient
convectionapplied
coefficient to the first to
applied 10thelayers
first of
10 layers o
the wall was assumed to be a constant value of 3 W/m 2 K corresponding to the estimated
the wall was assumed to be a constant value of 3 W/m K corresponding to the estimated
2

value for a vertical

valuesurface [63]. The
for a vertical ambient
surface [63].temperature
The ambientdefinition
temperature for definition
the first layer of the
for the first layer o
wall was set tothe74.5 ◦ C to approximate the air temperature near the heated bed, while the
wall was set to 74.5 °C to approximate the air temperature near the heated bed, while
ambient temperature
the ambientfor the remainingfor
temperature layers
the remaining 40 ◦ C,was
was set tolayers which
set tois40
approximately the
°C, which is approximately
value measured theby
value measured
embedded by embedded after
thermocouples thermocouples
cooling toaftersteady cooling
[Link] state. A convec
convection
coefficient of −tion
3 W/m 2 K was
coefficient ofimposed
−3 W/m2K at was imposedsurface
the bottom at the bottom
of the wallsurface of the wall
to obtain a “net to obtain a
convection” of “net
zeroconvection”
in that region. of zero in thatitregion.
Although Although itvalue,
is an unphysical is an this
unphysical
approach value, this approach
ensures
no convective ensures
behavior nobetween
convective thebehavior
part and between
bed. Thisthe is
part and [Link]
necessary Thistois limitations
necessary due of to limita
the Abaqus AM module, which is incapable of natively differentiating between exterior between
tions of the Abaqus AM module, which is incapable of natively differentiating
exteriorconvection
part surfaces where part surfaces where
occurs convection
and exterior occurs
faces of andtheexterior
part thatfaces areofinthe part that are in
contact
with the print bed. This negative convective coefficient approach could not be repeated fornot be re
contact with the print bed. This negative convective coefficient approach could
peated for condition
the radiation boundary the radiation thatboundary condition
is similarly [Link] isAssimilarly
a result,imposed.
the FE model As a result,
is the FE
model is expected to over‐predict cooling
expected to over-predict cooling at the bottom surface of the part. at the bottom surface of the part.
Thermal conductance values of 1, 5, 10, 25, 50, 75, and2 100 W/m2K were investigated
Thermal conductance values of 1, 5, 10, 25, 50, 75, and 100 W/m K were investigated.
Temperature data were extracted from two nodes at the interface over a period of 4000
Temperature data were extracted from two nodes at the interface over a period of 4000 s
seconds and compared against data reported by thermocouples C0 and C1. As shown in
and compared against data reported by thermocouples C02 and C1. As shown in Figure 5, a
Figure 5, a conductance (C) value of 10 W/m K produces thermal histories with minima
conductance (C) value of 10 W/m2 K produces thermal histories with minimal error. The
error. The difference of 5–10 °C at the steady‐state regime between experimental and
difference of 5–10 ◦ C at the steady-state regime between experimental and model data is
model data is deemed negligible in favor of the early period of cooling because it captures
deemed negligible in favor
the thermal of the early
behavior period
of layers closerof to
cooling
the [Link]
For both it captures
thermocouplesthe thermal
at the interface
behavior of layers closer to the bed. For both thermocouples at the interface,
temperature curves generated by using C = 10 W/m K exhibited the minimum tempera
2 temperature
curves generated 2 K exhibited the minimum temperature error
ture by
error using
untilC 750= seconds,
10 W/mcorresponding to approximately 5 layers of deposition in the
until 750 s, corresponding
manufacturing toprocess.
approximately
Figure 5 5shows
layersroot‐mean‐square
of deposition inerror the manufacturing
(RMSE) plots with a five
process. Figuresecond
5 shows root-mean-square error (RMSE) plots with
window, visualizing that C = 10 W/m2K is the conductance value a five-second window,
that minimizes
2
visualizing that C = 10 W/m K is the conductance value that minimizes error compared to
both thermocouples. Accordingly, the convective coefficient study utilizes C = 10 W/m2 K
as the thermal conductance value between the wall and the bed.
 Materials 2023, 16, x FOR PEER REVIEW 9 of 19

Materials 2023, 16, 6486 9 of 18

error compared to both thermocouples. Accordingly, the convective coefficient study uti‐
lizes C = 10 W/m2K as the thermal conductance value between the wall and the bed.

Figure 5. Temperature data comparison (top) and RMSE plots (bottom) between TCs C0 and C1
Figure 5. Temperature data comparison (top) and RMSE plots (bottom) between TCs C0 and C1 and
and FE model. The plots for C = 75 W/m2K and C = 100 W/m2K are indistinguishable.
FE model. The plots for C = 75 W/m2 K and C = 100 W/m2 K are indistinguishable.

[Link]. Wall Convection Study

Wall Convection Study
One of the most basic problems in the study of heat transfer over external surfaces is
One of the most basic problems in the study of heat transfer over external surfaces
the natural convection boundary layer flow over a semi‐infinite flat plate [65]. The free
is the natural convection boundary layer flow over a semi-infinite flat plate [65]. The
convection problem of a non‐isothermal vertical plate has been extensively studied by
free convection problem of a non-isothermal vertical plate has been extensively studied
several authors [66–69]. Analytical equations have been derived to obtain convection co‐
by several authors [66–69]. Analytical equations have been derived to obtain convection
efficients for a vertical plate through Nusselt number equations, which are a function of
coefficients for aand
the Rayleigh vertical
Prandtlplate through
numbers NusseltAlthough
[16,70,71]. number equations,
there is a small which airare a function
inlet near the of
theBAAM
Rayleigh
bed and Prandtl
to help remove numbers
vapors, [16,70,71].
the environmentAlthough
insidethere is a small
the printer air inlet
is assumed tonear
favorthe
BAAM bed to help remove vapors, the
free convection for the purposes of this study. environment inside the printer is assumed to favor
free convection for the purposes of this study.
The conductance between the wall and bed was set to 10 W/m K following the results
2
2 K following the results
of The conductance
the previous [Link] the wall
The ambient and bed was
temperature of set
the to 10 layer
first W/mof the wall was main‐
of the previous
tained at 74.5 study.
°C, whileThethe
ambient
ambienttemperature
temperature offor
thethe
first layer of the
remaining layerswall was
was setmaintained
to 40 °C
at in
74.5 ◦ C, while the ambient temperature for the remaining layers was set to240 ◦ C in
similar fashion to the conductance study. A convection coefficient (h) of −3 W/m K was
similar
againfashion
adoptedtofor the
theconductance
bottom surfacestudy. A wall.
of the convection
Severalcoefficient
iterations were −3 W/m2with
(h) ofcompleted K was
again
various convection coefficients of 3, 6, 9, 12, 15, 18, and 21 W/m K for the entire wall with
adopted for the bottom surface of the wall. Several iterations
2 were completed to
various convection
determine coefficients
the coefficient of 3, 6,
value which 9, 12, 15,the
minimized 18,error
andbetween
21 W/m 2 K for the
measured and entire
model‐wall
to predicted
determine temperatures.
the coefficientDatavalue
were which
extracted at the locations
minimized of all between
the error embeddedmeasuredthermocou‐ and
ples other than C3,
model-predicted AI0 and AI3, Data
temperatures. whichwere were extracted
excluded due at theto poor bonding.
locations of The temper‐
all embedded
ature distribution
thermocouples otherwas
thanextracted at nodes
C3, AI0 and AI3, that correspond
which were excluded to the due
physical location
to poor of theThe
bonding.
thermocouples
temperature and waswas
distribution recorded overatthe
extracted duration
nodes of the simulation.
that correspond Model‐predicted
to the physical location of
thetemperatures
thermocouples for and
various
wasconvection
recorded overcoefficient values are
the duration of theshown in FigureModel-predicted
simulation. 6.
Table 4for
temperatures shows that convection
various as the vertical position values
coefficient of eacharethermocouple
shown in Figure increases,
6. the con‐
vection coefficient that matched
Table 4 shows that as the vertical position of each thermocouple increases, These
experimental data best also tended to increase. re‐
the convec-
sults are corroborated via RMSE analysis. Accordingly, a coefficient
tion coefficient that matched experimental data best also tended to increase. These results value that minimized
areerror at each height
corroborated was selected
via RMSE [Link],
create a distribution of convection
a coefficient value thatcoefficients
minimizedthat error
minimized error over the entire wall. This progression of coefficients as a function of ver‐
at each height was selected to create a distribution of convection coefficients that minimized
tical distance from the print bed has also been observed with the use of IR cameras [15].
error over the entire wall. This progression of coefficients as a function of vertical distance
from the print bed has also been observed with the use of IR cameras [15].
Materials 2023, 16, 6486
Materials 2023, 16, x FOR PEER REVIEW 10 of 18 10 of 19

Figure 6.
Figure 6. Temperature Temperature
data comparisondata comparison
between between
TCs C2, TCsC6,
C4, C5, C2,C7,
C4,AI1,
C5, C6,
andC7,
AI2AI1,
and and
FEAI2 and FE model
model
with varying convection
with varying convection values. values.

Table 4. Best fitted value of convection coefficient h at each embedded thermocouple height for
Table 4. Best fitted value of convection coefficient (h) at each embedded thermocouple height for the
the FE model.
FE model.
Wall Height (m) Convection Coefficient—h (W/m2K)
Wall Height0.19
(m) Convection Coefficient—h
3 (W/m2 K)
0.19 0.39 3 6
0.39 0.58 6 9
0.58 0.78 9 12
0.78 0.98 12 15
0.98 15

3.3. Wall Residual Stress Study

3.3. Wall Residual Stress Study and stress distributions predicted by a finite element model assuming
Temperature
Temperature and stress
a constant distributions
convection predicted
coefficient by a finite
were compared element
against model
those assuming
utilizing a convec‐
a variable
constant convection coefficient
tion coefficient towere compared
determine againstonthose
the impact utilizing a of
the prediction variable
residualconvection
stresses generated
coefficient to determine
within the the impact
part. Modelon the prediction
geometry of residual
and material stresses
properties weregenerated
identical within
to those used in
the part. Model geometry and material properties were identical to those used in the
Materials 2023, 16, 6486 11 of 18

convective coefficient study. For the FE model using constant convection, the coefficient
was estimated analytically [70,71] by Equation (1).

Nu × k air
h= (1)
L
where kair is the thermal conductivity of the air, L is the wall height, and Nu is the average
Nusselt number. The average Nusselt number for a vertical isothermal plate is given by
Equation (2).
 2
 1/6 
 0.387 × Ra 
Nu = 0.825 + h i 8/27 
(2)
1 + (0.492/Pr )9/16

 

where Ra is the Rayleigh number (dimensionless number associated with buoyancy-driven

flow), Gr is the Grashof number (dimensionless number which approximates the ratio of
the buoyancy to viscous forces acting on a fluid), and Pr is the Prandtl number (ratio of
momentum diffusivity to thermal diffusivity), which are each defined by Equation (3),
Equation (4), and Equation (5), respectively.

Ra = Gr × Pr (3)

gβ(Tfilm − T∞ )L3
Gr = (4)
νair 2
νair
Pr = (5)
α air
In Equations (3)–(5), νair is the kinematic viscosity of the air, αair is the thermal dif-
fusivity of the air, g is the gravitational constant, β ≈ 1/T∞ is the volumetric expansion
coefficient of the air, Tfilm is the arithmetic mean between deposition temperature and T∞ ,
and T∞ is the ambient temperature. The properties of air used in Equations (3)–(5) [72] are
given in Table 5. The estimated constant convection coefficient value as calculated with
Equation (1) is 11.2 W/m2 K.

Table 5. Properties used in the calculation of the constant convection coefficient.

Properties Value
Ambient Temperature (◦ C) 40
Deposition Temperature (◦ C) 200
Wall Height (m) 1
Air Density (kg/m3 ) 1.127
Air Thermal Diffusivity (m2 /s) 2.346 × 10−5
Air Kinetic Viscosity (m2 /s) 1.702 × 10−5
Air Thermal Conductivity (W/mK) 0.02662
Gravitational Constant (m/s2 ) 9.81

Figure 7 shows a graph of the thermal comparison between FE models with constant
and variable convection coefficients. Temperatures were plotted for the nodes that coincide
with the dashed white line along the vertical axis of the wall at different moments during
the manufacturing process, as shown by images prepended to the top of the graph. Five
instances in time were captured corresponding to periods when the wall was at the same
height as each respective pair of embedded thermocouples. At Time 1 (1.7 h), the model
with constant convection over-predicts cooling by as much as 22 ◦ C, which is consistent
with the previous convection study showing that h = 3 W/m2 K is more suitable than
the area-averaged value of h = 11.2 W/m2 K. At Time 2 (3.6 h), the constant coefficient
model also over-predicts cooling, although the difference is reduced to 14 ◦ C at maximum.
Although at Time 3, there is still slight overcooling exhibited by the model with constant h,
Materials 2023,16,
Materials2023, 16, 6486
x FOR PEER REVIEW 12 of 1912 of 18

h,temperature
temperaturecurves
curvesatatTime
Time3 3and
andTime
Time4 4using
usingthe
theexperimentally
experimentallyderived,
derived,variable
variablevalues
values of h closely
of h closely approximate
approximate thosethose
basedbased
on theonconstant
the constant convection
convection coefficient
coefficient value.
value. At Time
At
5, Time 5, the with
the model modelconstant
with constant ℎ begins
h begins to under‐predict
to under-predict cooling,
cooling, whichwhich is expected
is expected because
because the experimentally
the experimentally derived derived coefficient
coefficient value value
for that that height
forheight 2 K) 2is
(15 W/m
(15 W/m K)higher
is higher
than the
than the constant
constant coefficient
coefficient value, value,
whichwhich
forcesforces heat dissipation
heat dissipation at a faster
at a faster [Link].

Figure
[Link] of temperature
Comparison distribution
of temperature across across
distribution the wallthe
height
wall(dashed
height white line)white
(dashed at dif‐line) at
ferent moments in time between the FE models with constant and variable convection coefficients.
different moments in time between the FE models with constant and variable convection coefficients.
The
Theresults
resultsfrom from both
boththermal
thermalmodels
modelswere given
were as inputs
given to structural
as inputs to structuralmodels for for
models
residual stress computation and comparison. Structural analyses were executed with an
residual stress computation and comparison. Structural analyses were executed with an
Abaqus AM module with the CTE and elastic stiffness constants determined through ma‐
Abaqus AM module with the CTE and elastic stiffness constants determined through mate-
terial characterization. The models utilized a 100 second time step and 8‐node linear hex‐
rial characterization. The models utilized a 100 s time step and 8-node linear hexahedral
ahedral elements with reduced integration with hourglass control (C3D8R). Mesh connec‐
elements with reduced integration with hourglass control (C3D8R). Mesh connectivity
tivity and nodal positions were retained from the thermal model. Fixed boundary condi‐
and nodal positions were retained from the thermal model. Fixed boundary conditions
tions were applied at the base of the wall to simulate a rigid adhesion to the bed. Stresses
were applied at
were measured the the
from base of the wall
integration to simulate
points a rigid
in elements adhesion
along to the
the height of thebed. Stresses
wall at the were
measured from the integration points in elements along the height
same instances in time that temperature data were extracted. Figure 8 shows the compar‐ of the wall at the same
instances
ison in time
of stresses that
in the temperature
vertical direction data were
(Szz) extracted.
between Figurewith
the models 8 shows
constantthe and
comparison
vari‐ of
stresses
able in the vertical
convection direction
coefficients. (Szz)component
This stress between the wasmodels
selected with constant
because and variable
it represents the con-
vectioninter‐layer
normal coefficients. This
stress stress beads,
between componentwhichwas selected
is a critical because when
parameter it represents
assessing thethe
normal
inter-layer
potential forstress between
inter‐layer beads, which
de‐bonding duringismanufacturing
a critical parameter
or afterwhen
a partassessing the potential
is completed.
for inter-layer de-bonding
A direct comparison during
between manufacturing
strains calculated by or the
after
FEamodel
part isandcompleted.
experimentally
A direct
obtained strains comparison
is not possible between
due tostrains calculated
the lack of availableby the
[Link]
such,and experimentally
the following
obtained
residual strains
stress is not possible
predictions are notdue to the lack of
experimentally available
validated. data. Asthe
However, such, the following
simulation‐
residual stress predictions are not experimentally validated.
based comparison shows that the model with constant h predicts higher magnitudes ofHowever, the simulation-based
Szz when compared
comparison shows that to thethevariable h model
model with with maximum
constant discrepancies
h predicts higher magnitudes in the of
lower
Szz when
regions
compared of theto wall. In Figure
the variable 8, the Szz
h model stress
with distributions
maximum show higher
discrepancies in thestress
lowerconcentra‐
regions of the
tion
[Link]
theFigure
base of 8,the
thewall
Szz for thedistributions
stress model with ashow constant convection
higher coefficient. at the base of
stress concentration
The PETG/CF
the wall for the model tensile strength
with was measured
a constant convection incoefficient.
accordance with ASTM D638 [73]
with anTheaverage
PETG/CF value for Szz of
tensile 9 MPa. Peak
strength was stress valuesinfrom
measured both constant
accordance withand ASTMvariable
D638 [73]
hwith
models were below
an average valuethisfor
threshold,
Szz of 9 which is consistent
MPa. Peak with the
stress values observation
from both constantthat noandcrack
variable
formation
h models wereor de‐bonding
below thisoccurred
threshold,in the
whichwall. Nevertheless,
is consistent with the
thecomparatively
observation that higherno crack
stresses
formationpresented by the constant
or de-bonding h model
occurred could
in the impact
wall. predictions the
Nevertheless, of crack formation inhigher
comparatively
astresses
part. In presented
summary, by thermal history and stress analysis results show that
the constant h model could impact predictions of crack formation neglecting the in
a part. In summary, thermal history and stress analysis results show that neglecting the
variation of h could impose inaccuracies on the quality of predictions generated by the
FE model.
 Materials 2023, 16, x FOR PEER REVIEW 13 of 19

Materials 2023, 16, 6486 variation of h could impose inaccuracies on the quality of predictions generated by the FE 13 of 18
model.

Figure 8. Comparison
Figure of stress
8. Comparison (Szz) (Szz)
of stress distribution across the
distribution wall the
across height (dashed
wall heightwhite line) white
(dashed at dif‐ line) at different
ferent moments in time between the FE models with constant and variable convection coefficients.
moments in time between the FE models with constant and variable convection coefficients.
3.4.
[Link] to Other
Extension Vertical
to Other Walls Printed
Vertical Walls on BAAMon BAAM
Printed
In the convective study presented previously, the convection coefficients for the
In the convective study presented previously, the convection coefficients for the
PETG/CF wall were found by means of comparing model and experimental temperature
data along thewall
PETG/CF were
height of thefound
wall. by
Themeans of comparing
same study was conducted model
for anand experimental
ABS/CF wall tempera-
ture data along the height of the wall. The same study was conducted
manufactured with the BAAM consisting of the same geometry. Different manufacturing for an ABS/CF wall
manufactured
parameters, such aswith
layerthe BAAM
time, consisting
were selected of the same
and accounted geometry.
for in the FE modelDifferent
throughmanufacturing
an updated event
parameters, series
such asinput.
layerThermocouple placementand
time, were selected in the ABS/CF wall
accounted forwas
in executed
the FE model through
inan
similar fashion to the PETC/CF equivalent. The objective for
updated event series input. Thermocouple placement in the ABS/CF repeating the study wall
in a was executed
different material system was to discover whether the ABS/CF case analysis would result
in similar fashion to the PETC/CF equivalent. The objective for repeating the study in a
in a similar h trend found for the PETG/CF wall, thus adding reliability to those results.
differentofmaterial
Repetition the studysystem
produced was to discover
similar whether
convection thevariations
coefficient ABS/CFascase analysis
presented in would result
a similar h
Table 6. Although the trend found for the ABS/CF wall study is not linear, the h values areto those results.
in trend found for the PETG/CF wall, thus adding reliability
Repetition
close of thetostudy
in magnitude produced
those reported similar
in Table convection
4 even when takingcoefficient variations
into consideration theas presented in
differences in material the
Table 6. Although systems
trendand printing
found forprocess conditions.
the ABS/CF wall study is not linear, the h values are
close in magnitude to those reported in Table 4 even when taking into consideration the
Table 6. Best fitted value of convection coefficient (h) for each inter‐layer TC height (h versus wall
differences
height) for the FEinmodel
material systems and printing process conditions.
of ABS/CF.

Wall Height (m) Convection Coefficient—h (W/m K) 2

Table 6. Best fitted value of convection coefficient (h) for each inter-layer TC height (h versus wall
0.19 3
height) for the FE model of ABS/CF.
0.39 3
0.58 Wall Height (m) 6 Convection Coefficient—h (W/m2 K)
0.78 9
0.98 0.19 15 3
0.39 3
0.58 6
The means of providing a generalized solution that could potentially be applied to
any vertical wall printed 0.78
on the BAAM is now presented. Researchers commonly 9 report
0.98 15

The means of providing a generalized solution that could potentially be applied to

any vertical wall printed on the BAAM is now presented. Researchers commonly report
generalized solutions through correlation equations, which are expressions developed
to provide approximations to local convection coefficients for specific conditions such
as vertical and horizontal plates, enclosures, etc. These equations are obtained through
numerical, analytical, and experimental measurement [70] methods. The assumption in
using correlation equations for natural convection in process simulations is that natural
convection occurs instantaneously during the manufacturing process [44]. Correlation
Materials 2023, 16, 6486 14 of 18

equations are generally posed in terms of dimensionless numbers such as the Rayleigh and
Prandtl numbers given by Equations (3) and (5).
Analytical solutions obtained for the laminar boundary layer problem of both isother-
mal and non-isothermal vertical plates provide expressions for the local Nusselt number
(Nuz ) in the same general format as Equation (6) [44,74]. In this study, the Nusselt number
was found by applying the relationship between (Nuz ) and h given by Equation (1), while
the Rayleigh number was found by using the surface temperature data as a function of
height from the FE model with the progression of convection coefficients that minimized
error in relation to experimental data from the PETG/CF wall. In Equation (4), before
applying Equation (3), the Prandtl number was simply calculated by Equation (5) given
the air properties according to the surface temperature data from the variable convection
coefficient wall model.
Nuz = A· Ram n
z · Pr z (6)
The unknown parameters A (521.22), m (0.26), and n (5.00) were found by utilizing
the nonlinear least-squares curve-fitting method in MATLAB. These parameters were then
used to derive convection coefficient values which were applied along the height of the
ABS/CF wall via Equation (1). This new distribution of h was then used in the FE model
of the ABS/CF wall, and the temperature results were compared with the equivalent
experimentally obtained thermocouple data. The comparison between thermocouple and
FE data shown in Figure 9 displays good agreement as corroborated by RMS analysis. The
result of this analysis demonstrates that the correlation equation derived from PETG/CF
data
Materials 2023, 16, x canREVIEW
FOR PEER potentially be applied to other vertical walls printed on the BAAM,15thus
of 19 extending

the solution found to situations having different material systems and process conditions.

Figure 9. Temperature
Figure 9. Temperature data comparison
data comparison betweenTCs
between TCs C2,
C2,C3,C3,
C4, C4,
C5, C6,
C5,C7,C6,
AI0,C7,
AI1,AI0,
AI2, and AI3AI2, and AI3
AI1,
and FE model using the variable convection coefficients found by Equation (6).
and FE model using the variable convection coefficients found by Equation (6).
4. Conclusions
This work successfully demonstrates how to obtain in situ temperature data from a
additively manufactured PETG/CF vertical wall using thermocouples embedded between
layers in the part. Combined with comprehensive material characterization, thermal his‐
tory data were given as inputs for the development of a predictive FE process model in
Abaqus CAE. The following can be concluded:
Materials 2023, 16, 6486 15 of 18

4. Conclusions
This work successfully demonstrates how to obtain in situ temperature data from a
additively manufactured PETG/CF vertical wall using thermocouples embedded between
layers in the part. Combined with comprehensive material characterization, thermal history
data were given as inputs for the development of a predictive FE process model in Abaqus
CAE. The following can be concluded:
1. There were significant differences in temperature predictions between FE models
that applied constant or variable convection coefficients. During the early stages of
manufacturing, the model with a constant convection coefficient over-predicts cooling
and gradually transitions to an under-prediction of cooling rate.
2. The progression of convection coefficients which minimized error in comparison with
experimental data increased with the height of the wall. This may be due to the
difference in air flow at regions of the wall at increasing heights above the BAAM
print bed. The linearity of the increase of convection coefficient is a compelling result.
However, this linear behavior cannot be extended as a rule, which is demonstrated by
independent predictions of the thermal history from the ABS/CF wall manufactured
with similar process conditions.
3. Significant differences were found in stress predictions between both models. The
model that used a constant convection coefficient predicted normal stresses in the
vertical direction at locations closest to the bed and during all stages of manufacturing
that were nearly double those predicted with the variable convection coefficient model.
Although cracks were not observed in the part and both models predicted stresses
below the material strength, these results suggest that accurate thermal modeling is
crucial for the prediction of residual stresses and, consequently, part integrity during
and immediately after the manufacturing process.
4. Regarding material characterization, adding thermal conductivity data as a function
of temperature to the FE thermal model may further improve model predictions.
5. The analysis of data generated by the PETG/CF wall experiment enabled the deriva-
tion of a correlation equation between convection coefficient and air properties that
can be applied to other vertical surfaces manufactured with the BAAM.

Author Contributions: Conceptualization, F.R.P., M.I., W.G.D. and R.A.L.-A.; methodology, F.R.P.,
M.I., W.G.D. and R.A.L.-A.; software, F.R.P. and M.I.; validation, F.R.P. and M.I.; formal analysis,
F.R.P. and M.I.; investigation, F.R.P. and M.I.; resources, R.A.L.-A.; data curation, F.R.P., M.I. and L.S.;
writing—original draft preparation, F.R.P.; writing—review and editing, F.R.P., M.I., L.S., W.G.D. and
R.A.L.-A.; visualization, F.R.P. and R.A.L.-A.; supervision, W.G.D. and R.A.L.-A.; project administra-
tion, W.G.D. and R.A.L.-A.; funding acquisition, R.A.L.-A. All authors have read and agreed to the
published version of the manuscript.
Funding: This research was funded by the U.S. Army Engineer Research and Development Center
(ERDC) through the Accelerated Rapid Prototyping (ARP) Project 1—Process and Modeling, Contract
No. W912HZ21C0005, the Russell S. Bodwell University Distinguished Professorship, and the
Malcolm G. Long ‘32 Professorship in Civil Engineering.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.
Acknowledgments: The authors want to thank Wesley Bisson, Richard Lafreniere Jr., and Audrey
Laffely for helping with the machining of the samples for material testing. The authors would also
like to thank George Scarlet and Britt Helten for providing valuable insights on FE modeling.
Conflicts of Interest: The authors declare that there are no conflict of interest, financial or personal,
regarding the publication of this paper.
Materials 2023, 16, 6486 16 of 18

References
1. Post, B.K.; Chesser, P.C.; Lind, R.F.; Roschli, A.; Love, L.J.; Gaul, K.T.; Sallas, M.; Blue, F.; Wu, S. Using Big Area Additive
Manufacturing to Directly Manufacture a Boat Hull Mould. Virtual Phys. Prototyp. 2019, 14, 123–129. [CrossRef]
2. Hassen, A.A.; Springfield, R.; Lindahl, J.; Post, B.; Love, L.; Duty, C.; Vaidya, U.; Pipes, R.B.; Kunc, V. The Durability of Large-
Scale Additive Manufacturing Composite Molds. In Proceedings of the CAMX Conference Proceedings, Anaheim, CA, USA,
27–29 September 2016.
3. Al Jassmi, H.; Al Najjar, F.; Mourad, A.-H.I. Large-Scale 3D Printing: The Way Forward. IOP Conf. Ser. Mater. Sci. Eng. 2018,
324, 012088. [CrossRef]
4. Holshouser, C.; Newell, C.; Palas, S.; Love, L.J.; Kunc, V.; Lind, R.F.; Lloyd, P.D.; Rowe, J.C.; Blue, C.A.; Duty, C.E.; et al. Out of
Bounds Additive Manufacturing. Adv. Mater. Process. 2013, 171, 15–17. [CrossRef]
5. Love, L.J.; Kunc, V.; Rios, O.; Duty, C.E.; Elliott, A.M.; Post, B.K.; Smith, R.J.; Blue, C.A. The Importance of Carbon Fiber to
Polymer Additive Manufacturing. J. Mater. Res. 2014, 29, 1893–1898. [CrossRef]
6. Love, L.J. Utility of Big Area Additive Manufacturing (BAAM) for the Rapid Manufacture of Customized Electric Vehicles; Oak Ridge
National Lab. (ORNL): Oak Ridge, TN, USA, 2015.
7. Duty, C.E.; Kunc, V.; Compton, B.; Post, B.; Erdman, D.; Smith, R.; Lind, R.; Lloyd, P.; Love, L. Structure and Mechanical Behavior
of Big Area Additive Manufacturing (BAAM) Materials. Rapid Prototyp. J. 2017, 23, 181–189. [CrossRef]
8. Ferrini-Mundy, J.; Varahramyan, K. 2022 Research Report: R1 Global Impact—Local Relevance; University of Maine: Orono, ME,
USA, 2023.
9. Lee, J.-Y.; An, J.; Chua, C.K. Fundamentals and Applications of 3D Printing for Novel Materials. Appl. Mater. Today 2017,
7, 120–133. [CrossRef]
10. Shofner, M.L.; Lozano, K.; Rodríguez-Macías, F.J.; Barrera, E.V. Nanofiber-reinforced Polymers Prepared by Fused Deposition
Modeling. J. Appl. Polym. Sci. 2003, 89, 3081–3090. [CrossRef]
11. Ahn, S.; Montero, M.; Odell, D.; Roundy, S.; Wright, P.K. Anisotropic Material Properties of Fused Deposition Modeling ABS.
Rapid Prototyp. J. 2002, 8, 248–257. [CrossRef]
12. Hoskins, D.; Kunc, V.; Hassen, A.; Lindahl, J.; Duty, C. Characterizing Thermal Expansion of Large-Scale 3D Printed Parts.
In Proceedings of the SAMPE 2019, Charlotte, NC, USA, 11 April 2019.
13. Advani, S.G.; Tucker, C.L. The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites. J. Rheol. 1987,
31, 751–784. [CrossRef]
14. Compton, B.G.; Post, B.K.; Duty, C.E.; Love, L.; Kunc, V. Thermal Analysis of Additive Manufacturing of Large-Scale Thermoplas-
tic Polymer Composites. Addit. Manuf. 2017, 17, 77–86. [CrossRef]
15. Brenken, B.; Barocio, E.; Favaloro, A.; Kunc, V.; Pipes, R.B. Development and Validation of Extrusion Deposition Additive
Manufacturing Process Simulations. Addit. Manuf. 2019, 25, 218–226. [CrossRef]
16. Brenken, B. Extrusion Deposition Additive Manufacturing of Fiber Reinforced Semi-Crystalline Polymers; Purdue University: West
Lafayette, IN, USA, 2017.
17. Zhang, Y.; Chou, Y.K. Three-Dimensional Finite Element Analysis Simulations of the Fused Deposition Modelling Process. Proc.
Inst. Mech. Eng. Part B J. Eng. Manuf. 2006, 220, 1663–1671. [CrossRef]
18. Lieneke, T.; Denzer, V.; Adam, G.A.O.; Zimmer, D. Dimensional Tolerances for Additive Manufacturing: Experimental Investiga-
tion for Fused Deposition Modeling. Procedia CIRP 2016, 43, 286–291. [CrossRef]
19. Alaimo, G.; Marconi, S.; Costato, L.; Auricchio, F. Influence of Meso-Structure and Chemical Composition on FDM 3D-Printed
Parts. Compos. Part B Eng. 2017, 113, 371–380. [CrossRef]
20. Ziemian, C.; Sharma, M.; Ziemian, S. Anisotropic Mechanical Properties of ABS Parts Fabricated by Fused Deposition Modelling.
Mech. Eng. 2012, 23, 159–180.
21. Soleyman, E.; Aberoumand, M.; Soltanmohammadi, K.; Rahmatabadi, D.; Ghasemi, I.; Baniassadi, M.; Abrinia, K.; Baghani, M.
4D Printing of PET-G via FDM Including Tailormade Excess Third Shape. Manuf. Lett. 2022, 33, 1–4. [CrossRef]
22. Soleyman, E.; Rahmatabadi, D.; Soltanmohammadi, K.; Aberoumand, M.; Ghasemi, I.; Abrinia, K.; Baniassadi, M.; Wang, K.;
Baghani, M. Shape Memory Performance of PETG 4D Printed Parts under Compression in Cold, Warm, and Hot Programming.
Smart Mater. Struct. 2022, 31, 085002. [CrossRef]
23. Dinwiddie, R.B.; Love, L.J.; Rowe, J.C. Real-Time Process Monitoring and Temperature Mapping of a 3D Polymer Printing. In
Proceedings of the SPIE 8705, Thermosense: Thermal Infrared Applications XXXV, Baltimore, MD, USA, 22 May 2013; Stockton, G.R.,
Colbert, F.P., Eds.; SPIE: Bellingham, WA, USA, 2013; Volume 8705, pp. 165–173.
24. Dinwiddie, R.B.; Kunc, V.; Lindal, J.M.; Post, B.; Smith, R.J.; Love, L.; Duty, C.E. Infrared Imaging of the Polymer 3D-Printing
Process. In Proceedings of the SPIE 9105, Thermosense: Thermal Infrared Applications XXXVI, Baltimore, MD, USA, 12 June 2014;
Colbert, F.P., Hsieh, S.-J., Eds.; SPIE: Bellingham, WA, USA, 2014; p. 910502.
25. Meraz Trejo, E.; Jimenez, X.; Billah, K.M.M.; Seppala, J.; Wicker, R.; Espalin, D. Compressive Deformation Analysis of Large Area
Pellet-Fed Material Extrusion 3D Printed Parts in Relation to in Situ Thermal Imaging. Addit. Manuf. 2020, 33, 101099. [CrossRef]
26. Planinsic, G. Infrared Thermal Imaging: Fundamentals, Research and Applications. Eur. J. Phys. 2011, 32, 1431. [CrossRef]
27. Kousiatza, C.; Tzetzis, D.; Karalekas, D. In-Situ Characterization of 3D Printed Continuous Fiber Reinforced Composites:
A Methodological Study Using Fiber Bragg Grating Sensors. Compos. Sci. Technol. 2019, 174, 134–141. [CrossRef]
Materials 2023, 16, 6486 17 of 18

28. Kousiatza, C.; Chatzidai, N.; Karalekas, D. Temperature Mapping of 3D Printed Polymer Plates: Experimental and Numerical
Study. Sensors 2017, 17, 456. [CrossRef] [PubMed]
29. Chin, R.K.; Beuth, J.L.; Amon, C.H. Thermomechanical Modeling of Molten Metal Droplet Solidification Applied to Layered
Manufacturing. Mech. Mater. 1996, 24, 257–271. [CrossRef]
30. Talagani, M.; DorMohammadi, S.; Dutton, R.; Godines, C.; Kumar Baid, H.; Abdi, F.; Kunc, V.; Compton, B.G.; Simunovic, S.; Duty,
C.; et al. Numerical Simulation of Big Area Additive Manufacturing (3D Printing) of a Full Size Car. SAMPE J. 2015, 51, 27–36.
31. Bhandari, S.; Lopez-Anido, R. Finite Element Analysis of Thermoplastic Polymer Extrusion 3D Printed Material for Mechanical
Property Prediction. Addit. Manuf. 2018, 22, 187–196. [CrossRef]
32. Bhandari, S.; Lopez-Anido, R.A. Discrete-Event Simulation Thermal Model for Extrusion-Based Additive Manufacturing of PLA
and ABS. Materials 2020, 13, 4985. [CrossRef] [PubMed]
33. Zhang, J.; Wang, X.Z.; Yu, W.W.; Deng, Y.H. Numerical Investigation of the Influence of Process Conditions on the Temperature
Variation in Fused Deposition Modeling. Mater. Des. 2017, 130, 59–68. [CrossRef]
34. Stockman, T.; Schneider, J.A.; Walker, B.; Carpenter, J.S. A 3D Finite Difference Thermal Model Tailored for Additive Manufactur-
ing. JOM 2019, 71, 1117–1126. [CrossRef]
35. Bhandari, S.; Lopez-Anido, R.A. Coupled Thermo-Mechanical Numerical Model to Minimize Risk in Large-Format Additive
Manufacturing of Thermoplastic Composite Designs. Prog. Addit. Manuf. 2023, 8, 393–407. [CrossRef]
36. Zhou, Y.; Nyberg, T.; Xiong, G.; Liu, D. Temperature Analysis in the Fused Deposition Modeling Process. In Proceedings of the
2016 3rd International Conference on Information Science and Control Engineering (ICISCE), Beijing, China, 8–10 July 2016;
pp. 678–682.
37. Ji, L.B.; Zhou, T.R. Finite Element Simulation of Temperature Field in Fused Deposition Modeling. Adv. Mater. Res. 2010, 97–101,
2585–2588. [CrossRef]
38. D’Amico, A.; Peterson, A.M. An Adaptable FEA Simulation of Material Extrusion Additive Manufacturing Heat Transfer in 3D.
Addit. Manuf. 2018, 21, 422–430. [CrossRef]
39. Bhandari, S.; Lopez-Anido, R. Finite Element Modeling of 3D-Printed Part with Cellular Internal Structure Using Homogenized
Properties. Prog. Addit. Manuf. 2019, 4, 143–154. [CrossRef]
40. Deering, R.A. Additive Manufacturing Part Level Distortion Sensitivity Analysis within Abaqus on a Thin Walled, Tubular
Structure. 2018. Available online: [Link]
PDF/2018-SaoE-Additive_Manufacturing_Part_Level_Distortion_Sensitivity_Analysis_within_Abaqus_on_a_Thin-walled_
_Tubular_Structure.pdf (accessed on 28 September 2023).
41. Cattenone, A.; Morganti, S.; Alaimo, G.; Auricchio, F. Finite Element Analysis of Additive Manufacturing Based on Fused
Deposition Modeling: Distortions Prediction and Comparison with Experimental Data. J. Manuf. Sci. Eng. 2019, 141, 011010.
[CrossRef]
42. Courter, B.; Savane, V.; Bi, J.; Dev, S.; Hansen, C.J. Finite Element Simulation of the Fused Deposition Modelling Process.
In Proceedings of the NAFEMS World Congress, Stockholm, Sweden, 11–13 June 2017; pp. 11–14.
43. Favaloro, A.J.; Brenken, B.; Barocio, E.; Pipes, R.B. Simulation of Polymeric Composites Additive Manufacturing Using Abaqus.
In Proceedings of the Science in the Age of Experience, Chicago, IL, USA, 15–18 May 2017; pp. 103–114.
44. Vaca, E.B. Fusion Bonding of Fiber Reinforced Semi-Crystalline Polymers in Extrusion Deposition Additive Manufacturing; Purdue
University: West Lafayette, IN, USA, 2018.
45. Nycz, A.; Lee, Y.; Noakes, M.; Ankit, D.; Masuo, C.; Simunovic, S.; Bunn, J.; Love, L.; Oancea, V.; Payzant, A.; et al. Effective
Residual Stress Prediction Validated with Neutron Diffraction Method for Metal Large-Scale Additive Manufacturing. Mater. Des.
2021, 205, 109751. [CrossRef]
46. Guessasma, S.; Belhabib, S.; Nouri, H. Printability and Tensile Performance of 3D Printed Polyethylene Terephthalate Glycol
Using Fused Deposition Modelling. Polymers 2019, 11, 1220. [CrossRef]
47. Dolzyk, G.; Jung, S. Tensile and Fatigue Analysis of 3D-Printed Polyethylene Terephthalate Glycol. J. Fail. Anal. Prev. 2019, 19,
511–518. [CrossRef]
48. Billah, K.M.M.; Lorenzana, F.A.R.; Martinez, N.L.; Wicker, R.B.; Espalin, D. Thermomechanical Characterization of Short Carbon
Fiber and Short Glass Fiber-Reinforced ABS Used in Large Format Additive Manufacturing. Addit. Manuf. 2020, 35, 101299.
[CrossRef]
49. Rios, O.; Carter, W.; Post, B.; Lloyd, P.; Fenn, D.; Kutchko, C.; Rock, R.; Olson, K.; Compton, B. 3D Printing via Ambient Reactive
Extrusion. Mater. Today Commun. 2018, 15, 333–336. [CrossRef]
50. Huang, H.; Ma, N.; Chen, J.; Feng, Z.; Murakawa, H. Toward Large-Scale Simulation of Residual Stress and Distortion in Wire
and Arc Additive Manufacturing. Addit. Manuf. 2020, 34, 101248. [CrossRef]
51. Nycz, A.; Kishore, V.; Lindahl, J.; Duty, C.; Carnal, C.; Kunc, V. Controlling Substrate Temperature with Infrared Heating to
Improve Mechanical Properties of Large-Scale Printed Parts. Addit. Manuf. 2020, 33, 101068. [CrossRef]
52. Colón Quintana, J.L.; Slattery, L.; Pinkham, J.; Keaton, J.; Lopez-Anido, R.A.; Sharp, K. Effects of Fiber Orientation on the Coeffi-
cient of Thermal Expansion of Fiber-Filled Polymer Systems in Large Format Polymer Extrusion-Based Additive Manufacturing.
Materials 2022, 15, 2764. [CrossRef] [PubMed]
53. ASTM D792–20; Standard Test Methods for Density and Specific Gravity (Relative Density) of Plastics by Displacement. West
ASTM: Conshohocken, PA, USA, 2020.
Materials 2023, 16, 6486 18 of 18

54. ASTM D3418–21; Standard Test Method for Transition Temperatures and Enthalpies of Fusion and Crystallization of Polymers by
Differential Scanning Calorimetry. ASTM: Conshohocken, PA, USA, 2021.
55. ISO/DIS 22007-2; Plastics—Determination of Thermal Conductivity and Thermal Diffusivity—Part 2: Transient Plane Source
Method. ISO: Geneva, Switzerland, 2022.
56. ASTM E831-19; Standard Test Method for Linear Thermal Expansion of Solid Materials by Thermomechanical Analysis. ASTM:
Conshohocken, PA, USA, 2019.
57. ASTM D7028-07; Standard Test Method for Glass Transition Temperature (DMA Tg) of Polymer Matrix Composites by Dy-Namic
Mechanical Analysis (DMA). ASTM: Conshohocken, PA, USA, 2015.
58. ASTM D5023-15; Standard Test Method for Plastics: Dynamic Mechanical Properties: In Flexure (Three-Point Bending). ASTM:
Conshohocken, PA, USA, 2016.
59. Steva, B.; Warren, K.; Seigars, C.; Helten, B. Evaluation of the Shear Behavior of CF-PETG and CF-PC Coupons Manufactured
Using Large-Scale Additive Manufacturing Processes. In Proceedings of the ECCM20—The 20th European Conference on
Composite Materials, Lausanne, Switzerland, 26–30 June 2022.
60. Seigars, C.; Warren, K.; Steva, B.; Murphy, C.; Helten, B. Characterizing the Tensile and Compressive Behavior of PETG/CF and
PC/CF Manufactured Using Large Scale Additive Processes. In Proceedings of the ECCM20—The 20th European Conference on
Composite Materials, Lausanne, Switzerland, 26–30 June 2022.
61. Kim, P.; Baid, H.; Hassen, A.; Kumar, A.; Lindahl, J.; Hoskins, D.; Ajinjeru, C.; Duty, C.; Yeole, P.; Vaidya, U.; et al. Analysis on
Part Distortion and Residual Stress in Big Area Additive Manufacturing with Carbon Fiber-Reinforced Thermoplastic Using
Dehomogenization Technique. In Proceedings of the Composites and Advanced Materials Expo (CAMX 2019), Anaheim, CA,
USA, 23–26 September 2019.
62. Choo, K.; Friedrich, B.; Daugherty, T.; Schmidt, A.; Patterson, C.; Abraham, M.A.; Conner, B.; Rogers, K.; Cortes, P.; MacDonald, E.
Heat Retention Modeling of Large Area Additive Manufacturing. Addit. Manuf. 2019, 28, 325–332. [CrossRef]
63. Morgan, R.V.; Stowers Reid, R.; Baker, A.M.; Lucero, B.; Bernardin, J.D. Emissivity Measurements of Additively Manufactured
Materials; Los Alamos National Lab. (LANL): Los Alamos, NM, USA, 2017.
64. Na, T.Y. Numerical Solution of Natural Convection Flow Past a Non-Isothermal Vertical Flat Plate. Appl. Sci. Res. 1977, 33,
519–543. [CrossRef]
65. Chen, T.Y.W.; Wollersheim, D.E. Free Convection at a Vertical Plate with Uniform Flux Condition in Non-Newtonian Power-Law
Fluids. J. Heat Transf. 1973, 95, 123–124. [CrossRef]
66. Sparrow, E.M.; Gregg, J.L. Laminar Free Convection from a Vertical Plate with Uniform Surface Heat Flux. J. Fluids Eng. 1956, 78,
435–440. [CrossRef]
67. Touloukian, Y.S.; Hawkins, G.A.; Jakob, M. Heat Transfer by Free Convection from Heated Vertical Surfaces to Liquids. J. Fluids
Eng. 1948, 70, 13–17. [CrossRef]
68. Cheesewright, R. Turbulent Natural Convection from a Vertical Plane Surface. J. Heat Transf. 1968, 90, 1–6. [CrossRef]
69. Churchill, S.W.; Chu, H.H.S. Correlating Equations for Laminar and Turbulent Free Convection from a Vertical Plate. Int. J. Heat
Mass Transf. 1975, 18, 1323–1329. [CrossRef]
70. Incropera, F.P.; DeWitt, D.P.; Bergman, T.L.; Lavine, A.S. Fundamentals of Heat and Mass Transfer, 6th ed; Wiley: New York, NY,
USA, 2007.
71. Çengel, Y.A.; Ghajar, A.J. Heat and Mass Transfer: Fundamentals and Applications; McGraw-Hill Education: New York, NY, USA, 2020.
72. Properties of Air at 1 Atm Pressure. Engineers Edge. Available online: [Link]
air_at_1_atm_pressure_13828.htm (accessed on 28 September 2023).
73. ASTM D638-22; Standard Test Method for Tensile Properties of Plastics. ASTM: Conshohocken, PA, USA, 2022.
74. Havet, M.; Blay, D. Natural Convection over a Non-Isothermal Vertical Plate. Int. J. Heat Mass Transf. 1999, 42, 3103–3112.
[CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.