Abstract: The flow drag of a particle in unbounded uniform incoming flow is one of the fundamental problems in particulate multi-phase flow. Particles in relevant flows are generally non-spherical. Earlier investigations used the drag of a volume-equivalent sphere as a first approximation. Recently more empirical formulas, mostly derived from data fitting of measurements, have been developed. In this paper, high precision numerical simulations are carried out for a fixed oblate/prolate spheroid in a flow with varying particle Reynolds number from 0.1 to 1000, where different spheroid orientation has been taken into account. The simulation model and method have been verified against known analytical solutions. For more general flow conditions without known analytical solutions, numerical results have been compared with those from various popular empirical formulas. Generally speaking, a spheroid experiences the highest drag when it is located in the flow with the broadest side against the flow. However, numerical simulation at Reynolds number 1000 reveals an exception, i.e., an oblate spheroid experiences higher drag when put slanted in the flow than when put crosswise. This is because the flow around the spheroid is more stable in an axial-symmetrical setting than in the asymmetrical setting when the wake becomes chaotic and hence the form drag increases substantially. Overall, the numerical simulation is believed more precise in predicting the drag of a particle of specific shape than empirical formulas which aim to apply for a broad range of shapes. Through accumulating more high precision numerical data, it is promising to develop a more precise and broadly applicable drag formula.