In this work, we present a complete methodology aiming to simulate the industrial process of polyurethane foams injection – expansion. Our method is based on both the use of strong numerical tools as level-set or auto adaptive anisotropic meshing (AAA), and a realistic description of all physical phenomena that occur. To the classical Navier-Stokes and heat equations are added a reaction equation for polymerization and a reaction equation for foaming. All these equations are coupled together as the heat fasten foaming reaction and polymerization, foaming reaction generates volumetric expansion, polymerization leads to viscosity freezing and slows down foaming, … The originality of our approach lies in enhancing each equation stabilization with a global stabilization scheme, ensuring an optimized coherency between these strongly coupled equations. Furthermore, behavior of liquid polyurethane during injection step has a great impact on final foam distribution. The use of multicriteria AAA meshing with a level-set method enables us to simulate the flow inside the cavity. This coupled with the stabilized governing equations enable to compute in a single step injection and expansion, even for large industrial parts. All these method used together enable us to give a realistic description of foam distribution and three dimensional density prediction. Finally, several analytical and industrial cases help us to validate the whole method.