## 2 Development of Beam Element-Based Simulation  
### Model Description  
A two-dimensional, beam-based, small-strain and deformation finite element model is presented here. As discussed in the introduction, an FDM part consists of layers made up of one or more beads, which all bond to their neighbors. As shown in Figure 1, the beads are represented in the model with beam elements and the inter-bead bond is simulated using multipoint constraints and vertical axial-only elements (termed vertical springs for the rest of this paper). The model is structural only and requires temperatures versus location and time as input. Examining the FDM process, the print head continuously deposits molten material (which forms a bead) at a given temperature, moving at a certain speed and extrusion rate, pauses to allow the bead to cool if directed by the G-code and then repeats the process for the next bead. This continuous process is approximated using discrete time steps in the model and in each time step only introducing elements whose centroids have already been printed. The settings defining this process determine the cross-section of the bead in addition to the area of the inter-bead bond, and the aspects relevant to the model are captured via the elemental cross-sectional areas (\(A_{bead}\) or \(A_{vert}\)) and moment of inertia (\(I_{bead}\) or \(I_{vert}\)). As the bead's temperature fluctuates due to cooling and rewarming due to the placement of adjacent beads, the stiffness of the bead increases and decreases, respectively. This phenomenon is incorporated into the model by defining the elastic modulus (\(E_{bead}\) or \(E_{vert}\)) as a function of temperature. Changes in bead temperature also induce thermal strains, which are included in the model as element pre-strains. It is assumed that thermal strains  
Figure 1: Element discretization. accumulate only once the bead temperature falls below the glass transition temperature of the polymer (\(T_{g}\)). The stiffness of the inter-bead bond is dependent on printing parameters and is accounted for in the vertical direction in the model by allowing \(E_{vert}\
eq E_{bead}\), where both moduli of elasticity can be calibrated to the specific printer and parameters through material testing. Linear, multipoint constraints are employed to enforce inter-bead shear transfer and assume full composite action (i.e., no interlayer slip). Flowcharts for both the process and the model's representation of that process are shown in Figure 2. The grey dotted lines denote correspondence between process and model steps. All aspects of the model are elaborated on in the following sections. Employing the principle of virtual work for the current model produces Equation (1),  
\[\delta W_{int_{b}}+\delta W_{int_{ax}}+\delta W_{int_{f}}-\delta W_{ext}=0 \tag{1}\]  
where sub-subscript \(b\) denotes bending, \(ax\) denotes axial and \(f\) denotes a spring foundation supporting the base of the printed part, which is discussed later. Substituting the well-established expressions for each virtual work term produces the weak form of the governing differential equation for a beam element in the model, shown in Equation (2),  
\[\int_{0}^{l}EIv^{\prime\prime}\delta v^{\prime\prime}dx+\int_{0}^{l}A\sigma \delta edx+\int_{0}^{l}kv\delta vdx-\int_{0}^{l}q\delta vdx=0 \tag{2}\]  
Figure 2: Process-model flowchart. where \(E\) is modulus of elasticity; \(I\) is the second area moment of inertia; a prime denotes differentiation with respect to \(x\), the position along the element; \(v\) is the transverse displacement along the element; \(\delta\) denotes the variation of a displacement, strain or curvature; \(A\) is the cross-sectional area of the element; \(\sigma\) is the axial stress; \(\epsilon\) is the axial strain; \(k\) is the stiffness of the foundation; \(q\) represents the applied vertical loads due to self-weight; \(u\) is the axial displacement and \(U^{el}\) is a vector of element nodal displacements. Thermal pre-strains are incorporated in the axial stress-strain relationship.