A smoothed particle hydrodynamics (SPH) numerical method is first developed to solve thermal transport problems in heterogeneous materials with finite thermal contact conductance and discontinuous temperature field at interfaces between individual grain-like or fiber-like particles. It is applied to study thermal transport and calculate effective conductivity in a linear chain of powder grains, two-dimensional, and three-dimensional random powder systems composed of spherical grains and high-aspect-ratio spherocylinders, and in a nanocomposite material with carbon nanotubes and polymer matrix. In all cases, the numerical errors are found to be sufficiently small at relatively large spacing between SPH particles. The SPH method is further developed to account for melting and solidification of the target material based on the enthalpy formulation and major interfacial effects, including surface tension force, Marangoni stresses, recoil force due to vapor pressure, mass removal, and evaporative cooling and is coupled with the ray tracing (RT) method to simulate the propagation, multiple reflections and absorption of incident laser radiation. The SPH-RT method is used to investigate the melt pool expulsion mechanism during high aspect ratio laser drilling of aluminum, and 316L stainless steel bulk targets. The drilling velocity is found to be non-linear which is dependent on the shape of the heating surface. The simulations reveal very strong effect of multiple reflections of laser radiation inside the keyhole on the drilling velocity. The main driving force of melt expulsion is the repulsive force produced by vapor pressure, while the Marangoni stresses only marginally affect the drilling velocity. The spattering of melt pool occurs at high laser intensities which is well captured by SPH-RT method. For small intensity a steady state is achieved where the incident laser energy is balanced by the thermal conduction and evaporative cooling and hence the drilling stops. The SPH-RT method is also used to study the drilling of copper with bursts of nanosecond laser pulses. The ablated material, surface temperature, and cavity depth after first pulse match with the solution obtained with the finite-difference solution of the heat conduction equation. The rate of ablation increases with number of laser pulses.