•Unified single-node curved boundary conditions based on TLBM are proposed with higher accuracy and parallel scalability than conventional schemes.•Parametric combined momentum-exchange method accurately simulates the lift coefficients for airfoils in transitional flows using Multi-GPU computation.•Simulated average Nussel numbers Nu¯ for forced convection of airfoils similar to the experimental results are equivalent to and modify the Hilpert correlation. Laminar separation bubbles in transitional flows of low-Reynolds-number airfoils significantly affect aerodynamic efficiency and heat transfer characteristics, which can be dealt with using the numerical method with boundary treatment of complex geometry. This work proposes the unified single-node curved boundary conditions based on the thermal lattice Boltzmann method, and several benchmarks verify the higher accuracy and parallel scalability over conventional schemes. Multi-GPU simulations for airfoil Eppler 61 and SD8020 are carried out, and the lift coefficients CL are accurately simulated by the parametric combined momentum-exchange method. Furthermore, time-averaged pressure coefficients Cp and local Nussel number Nu of airfoil surfaces are obtained, and the unsteady flow and heat transfer characteristics of laminar separation bubbles are captured. Simulated average Nussel numbers Nu¯ for an Eppler 61 airfoil with different Reynold numbers are equivalent to and modify the Hilpert correlation. The numerical method can quickly and accurately predict transitional flows of low-Reynolds-number airfoils and guide the design and application of heat transfer with blade configuration.