Addressing challenges in urban water infrastructure systems, including aging infrastructure, supply uncertainty, extreme events, and security threats, depends highly on water distribution networks modeling emphasizing the importance of realistic assumptions, modeling complexities, and scalable solutions. In this study, we propose a derivative‐free, linear approximation for solving the network water flow problem. The proposed approach takes advantage of the special form of the nonlinear head loss equations, and, after the transformation of variables and constraints, the water flow problem reduces to a linear optimization problem that can be efficiently solved by modern linear solvers. Ultimately, the proposed approach amounts to solving a series of linear optimization problems. We demonstrate the proposed approach through several case studies and show that the approach can model arbitrary network topologies and various types of valves and pumps, thus providing modeling flexibility. Under mild conditions, we show that the proposed linear approximation converges. We provide sensitivity analysis and discuss in detail the current limitations of our approach and suggest solutions to overcome these. All the codes, tested networks, and results are freely available on Github for research reproducibility.