Two nonlinear κ-ε models with the wall function method are applied to the fully developed turbulent flow in a square duct. Typical predicted quantities such as axial and secondary velocities, turbulent kinetic energy and Reynolds stresses are compared in details both qualitatively and quantitatively with each other. A nonlinear κ-ε model with the wall function method capable of predicting accurately duct flows involving turbulence-driven secondary motion is presented in the present paper. The nonlinear κ-ε model of Shih et al.[1] adopted in a commercial code is found to be unable to predict accurately duct flows with the prediction level of secondary flows one order less than that of the experiment.