Aircraft in-flight icing is a major safety issue for civil aviation, having already caused hundreds of accidents and incidents related to aerodynamic degradation due to post takeoff ice accretion. Airplane makers have to protect the airframe critical surfaces against ice build up in order to ensure continued safe flight. Ice protection is typically performed by mechanical, chemical, or thermal systems. One of the most traditional and still used techniques is the one known as hot-air anti-icing, which heats the interior of the affected surfaces with an array of small hot-air jets generated by a piccolo tube. In some cases, the thermal energy provided by hot-air ice protection systems is high enough to fully evaporate the impinging supercooled droplets (fully evaporative systems), while in other cases, it is only sufficient to maintain most of the protected region free of ice (running wet systems). In the latter case, runback ice formations are often observed downstream of the wing leading edge depending on hot-air, icing, and flight conditions. The design process of hot-air anti-icing systems is traditionally based on icing wind tunnel experiments, which can be very costly. The experimental effort can be significantly reduced with the use of accurate three-dimensional computational fluid dynamic (CFD) simulation tools. Nevertheless, such type of simulation requires extensive CPU time for exploring all the design variables. This thesis deals with the development of an efficient hot-air anti-icing system simulation tool that can reduce the computational time to identify the critical design parameters by at least two orders of magnitude, as compared to 3-d CFD tools, therefore narrowing down the use of more sophisticated tools to just a small subset of the entire design space. The hot-air anti-icing simulation tool is based on a combination of available CFD software and a thermodynamic model developed in the present work. The computation of the external flow properties is performed with FLUENT (in a 2-d domain) by assuming an isothermal condition to the airfoil external wall. The internal skin heat transfer is computed with the use of local Nusselt number correlations developed through calibrations with CFD data. The internal and external flow properties on the airfoil skin are provided as inputs to a steady state thermodynamic model, which is composed of a 2-d heat diffusion model and a 1-d uniform film model for the runback water flow. The performance of the numerical tool was tested against 3-d CFD simulation and experimental data obtained for a wing equipped with a representative piccolo tube anti-icing system. The results demonstrate that the simplifications do not affect significantly the fidelity of the predictions, suggesting that the numerical tool can be used to support parametric and optimization studies during the development of hot-air anti-icing systems.