Numerical simulations of multiphase flows have recently expanded beyond Newtonian fluids to include complex interfacial physics such as viscoelasticity and surfactant effects. Foam flows composed of many interacting bubbles are commonly encountered in industrial and environmental applications, including foam flooding for enhanced oil recovery in porous media. However, their numerical simulation remains highly challenging due to the presence of thin liquid films, strong interfacial deformation, and multi-bubble interactions. To accurately resolve gas–liquid interfaces and liquid films, we employ a high-resolution multiphase flow solver based on a weakly compressible scheme combined with Adaptive Mesh Refinement (AMR), which achieves good scalability on multi-GPU platforms[1]. To prevent unphysical bubble coalescence and collapse in dense foam systems, we adopt a Multi-Phase Field (MPF) approach in which each bubble is assigned an individual phase field, enabling stable simulations of multiple interacting bubbles. Furthermore, the stability of liquid films strongly depends on surfactant transport and interfacial adsorption/desorption processes. While conventional finite-volume-based surfactant transport models suffer from large numerical diffusion when realistic diffusion coefficients are used, we have developed a novel phase-field-based surfactant transport model[2] that preserves concentration while maintaining numerical stability. We demonstrate the capability of the proposed framework through simulations of surfactant-stabilized foam flows, highlighting the roles of interfacial transport, bubble interactions, and liquid film stability. The combination of AMR-based high-resolution simulations, the MPF method, and the new surfactant transport model provides a robust tool for investigating complex multiphase flow phenomena involving multiple deformable interfaces.