Face milling commonly generates surface quality of variation, is especially severe for milling of large-scale components with complex surface geometry such as cylinder block, engine head, and valve body. Thus surface variation serves as an important indicator both for machining parameter selection and components' service performance such as sealing, energy consumption, and emission. An efficient and comprehensive numerical model is highly desired for the prediction of surface variation of entire surface. This study proposes a coupled numerical simulation method, updating finite element (FE) model iteratively based on integration of data from abaqus and matlab, to predict surface variation induced by face milling of large-scale components with complex surfaces. Using the coupled model, three-dimensional (3D) variation of large-scale surface can be successfully simulated by considering face milling process including dynamic milling force, spiral curve of milling trajectory, and intermittently rotating contact characteristics. Surface variation is finally represented with point cloud from iterative FE analysis and verified by face milling experiment. Comparison between measured and predicted results shows that the new prediction method can simulate surface variation of complex components well. Based on the verified model, a set of analyses are conducted to evaluate the effects of local stiffness nonhomogenization and milling force variation on machined surface variation. It demonstrates that surface variation with surface peaks and concaves is strongly correlated with local stiffness nonhomogenization especially in feed direction. And thus the coupled prediction method provides a theoretical and efficient way to study surface variation induced by face milling of large-scale complex components.