Numerical models of heat conduction and natural convection are developed based on Meshless Moving particle semi-implicit (MPS) method to study the unsteady complex flow coupling with heat transfer. Fundamental heat exchange process is systematically investigated in one dimension and extended to two dimensions without internal heat source. The energy conservation equation is discretized on particles and coupled with flow by momentum change that caused by temperature difference. Both Dirichlet and Neumann boundary conditions are proposed with mirror image technique and the numerical results are verified to have high accuracy against benchmark solutions (theoretical or experimental, and or other numerical solutions) for several representative examples. To improve the calculation efficiency, symmetry boundary models are established and applied in transient heat conduction and natural convection heat transfer. The results using the symmetry boundary are comparable and reliable to the full-domain results, both of which agree well with the experimental and other numerical solutions. (C) 2015 Elsevier Ltd. All rights reserved.