Fused deposition modeling (FDM) is one of the most affordable and widespread additive manufacturing (AM) technologies. Despite its simplistic implementation, the physics behind this FDM process is very complex and involves rapid heating and cooling of the polymer feedstock. As a result, highly non-uniform internal stresses develop within the part, which can cause warpage deformation. The severity of the warpage is highly dependent on the process parameters involved, and therefore, currently extensive experimental studies are ongoing to assess their influence on the final accuracy of the part. In this study, a thermomechanical Finite Element model of the 3D printing process was developed using ANSYS. This model was compared against experimental results and several other analytical models available in the literature. The developed Finite Element Analysis (FEA) model demonstrated a good qualitative and quantitative correlation with the experimental results. An L9 orthogonal array, from Taguchi Design of Experiments, was used for the optimization of the warpage based on experimental results and numerical simulations. The optimum process parameters were identified for each objective and parts were printed using these process parameters. Both parts showed an approximately equal warpage value of 320 μm, which was the lowest among all 10 runs of the L9 array. Additionally, this model is extended to predict the warpage of FDM printed multi-material parts. The relative percentage error between the numerical and experimental warpage results for alternating and sandwich specimens are found to be 1.4% and 9.5%, respectively.