Few studies jointly investigate thermal and turbulent effects. In general, these subjects are treated separately. The purpose of this paper is to use the Immersed Boundary Method (IBM) coupled with the Virtual Physical Model (VPM) to investigate incompressible two-dimensional Newtonian flow around a heated square cylinder at constant temperature on its surface with forced convection and turbulence. The VPM model dynamically evaluates the force that the fluid exerts on the immersed surface and the thermal exchange between both in the Reynolds numbers (Re) window 40 ≤ Re ≤ 5×10³ . For simulations of turbulence the Smagorinsky and Spalart-Allmaras models are used. The first model uses the Large Eddy Simulation (LES) methodology and is based on the local equilibrium hypothesis for small scales associated with the Boussinesq hypothesis, such that the energy injected into the spectrum of the turbulence balances the energy dissipated by convective effects. The second model uses the concept Unsteady Reynolds Averaged Navier-Stokes Equations (URANS), with only one transport equation for turbulent viscosity, being calibrated in pressure gradient layers. The goal of this work is to analyse the combination of the heat-transfer phenomena with the turbulence for the thermo-fluid-structure interaction in a square cylinder. For this, it was developed a C/C++ code that requires low computational costs in regards to memory and computer facilities. It is observed that, with the increase of the Reynolds number, an increase of the drag coefficient occurs, as well as reinforces the influence of the pressure distribution downstream of the cylinder, which is strongly influenced by the formation and detachment of vortices on the upper and lower sides of the square cylinder.