A new numerical method is proposed to study two-phase flow and heat transfer for interlayer cooling of the new generation of multi-stacked computer chips. The fluid flow equations are developed in 3-dimensions based on the Arbitrary Lagrangian-Eulerian formulation (ALE) and the Finite Element Method (FEM), creating a new two-phase method with an improved model for the liquid-gas interface. A new adaptive mesh update procedure is also proposed for effective management of the mesh at the two-phase interface to remove, add and repair surface elements, since the computational mesh nodes move according to the flow. The Lagrangian description explicitly defines the two-phase interface position by a set of interconnected nodes which ensures a sharp representation of the boundary, including the role of the surface tension. The new methodology for computing the curvature leads to accurate results with moderate programming effort and computational cost. Static and dynamic tests have been carried out to validate the method and so far all the obtained results have compared well to analytical solutions and experimental results found in the literature, demonstrating that the new proposed methodology to simulate two-phase flows provides good accuracy to describe the interfacial forces and bubble dynamics. The new code was then used to simulate elengated bubble flows in square microchannels, being considered for two-phase interlayer cooling in future 3D-IC compute chips.