A theoretical model of thermal boundary layers and acoustic heating in microscale acoustofluidic devices is presented. Based on it, an iterative numerical model is developed that enables numerical simulation of nonlinear thermoviscous effects due to acoustic heating and thermal advection. Effective boundary conditions are derived and used to enable simulations in three dimensions. The theory shows how friction in the viscous boundary layers causes local heating of the acoustofluidic device. The resulting temperature field spawns thermoacoustic bulk streaming that dominates the traditional boundary-driven Rayleigh streaming at relatively high acoustic energy densities. The model enables simulations of microscale acoustofluidics with high acoustic energy densities and streaming velocities in a range beyond the reach of perturbation theory, and is relevant for design and fabrication of high-throughput acoustofluidic devices.