In this work we investigate the role of elasticity on the thermal Marangoni migration of a Newtonian droplet surrounded by a viscoelastic fluid matrix in a three‐ dimensional geometry for the case of small Reynolds and Marangoni numbers. The study has been conducted in the framework of a coupled Level‐Set‐Volume of Fluid method implemented using the CFD toolbox OpenFOAM. The resulting approach was validated in a variety of flow conditions by comparing our results with analytical correlations and relevant experimental data available in literature. In the present numerical experiments, we consider a neutrally buoyant system of a Newtonian droplet placed in a container with square cross‐section filled with an Oldroyd‐B fluid (a viscoelastic fluid of constant shear viscosity). We apply a thermal gradient by keeping two sides of the box at a different constant temperature so that the temperature gradients at the liquid‐liquid interface generate an imbalance in the interfacial stresses. Such imbalance in turn is responsible of the motion of the fluid from the higher temperature region to the lower temperature region. This mechanism results in the drop moving in the opposite direction due to the thrust generated by the counter motion of the surrounding phase. In order to quantify the viscoelastic effects, we introduce a new dimensionless parameter measuring the relative importance of thermocapillary and elastic stresses. According to the numerical results, the droplet migration speed and shape are significantly different from those observed for the Newtonian‐Newtonian system. This departure of the observed dynamics from Newtonian behaviour can be ascribed to the complex interplay between different effects, including droplet morphological evolution and related distribution of surface‐tension‐driven and elastic stresses at the interface
