We present a numerical homogenization procedure for nonlinear Darcy flow for laminar flow through arrays of generic, simply-connected, 2D microstructures. Machine learning is utilized to fit coefficients in the nonlinear law based on metrics computed from the solutions to a series of well-posed cell problems. Included in these coefficients is a mean flow-dependent nonlinear correction term which modifies the typical linear Darcy's law to account for inertial effects of the flow. The proposed method differs from existing works on fluid homogenization in several key ways. First, the nonlinear correction term enables accurate modeling of fluid behavior for Reynolds numbers up to, and even beyond, Re = 100. Second, the method is founded on a rigorously defined and derived nonlinear correction, rather than correlations or empirical models, the latter of which restrict the usefulness of the underlying scheme. Finally, our method is formulated to be sufficiently flexible so as to allow for the consideration of a wide variety of microstructures. Because of these features, our method can be applied to a range of practical analysis and design problems, yielding efficient approximations without the need for the full-scale resolution of individual microstructures. The correctness and efficacy of our numerical implementation is demonstrated via systematic comparisons with published experimental and numerical results, as well as full-scale analyses. In each case, it is shown that the nonlinear Darcy's Law model is able to achieve significantly greater accuracy than the typical linear Darcy's Law model as the Reynolds number of the flow increases beyond one, for little increase in computational cost.