The determination of the rotational speed and massflow of the fan of a turbofan at windmill is critical in the design of the engine-supporting structure and the sizing of the vertical stabilizer. Given the very high bypass ratio obtained at windmill, the flow in the fan stage and bypass duct is of prime interest. Classical CFD simulations have been shown to predict such flows accurately, but extensive parametric studies can be needed, stressing the need for reduced-cost modeling of the flow in the engine. A Body Force Modeling (BFM) approach for windmilling simulations is examined in the present contribution. The BFM approach replaces turbomachinery rows by source terms, reducing the computational cost (here by a factor 5). A shaft model is coupled to the BFM source terms, to drive the simulation to a power balance of the low-pressure shaft. The overall approach is thus self-contained and can predict both the massflow and the rotational speed in the windmilling regime. Comparisons with engine experimental results show the proposed model can predict the rotational speed within 7 %, and the massflow within 5 %. Local analysis and comparisons with experimental data and reference blade calculations show that the work exchange, in term of total temperature variation, is predicted within 0.5 K, and the overall total pressure ratio within 1 %. However, the losses in the stator are largely underestimated, which explains the discrepancy for the massflow predictions.