The analysis of dynamic soil-pile interaction problems requires the relation of soil resistance to lateral loading that is represented by nonlinear p-y curves in the beam on the nonlinear Winkler foundation (BNWF) approach. Current methods for p-y curves are either based on static load tests, or dynamic tests, which cannot properly consider soil nonlinearity. This study investigates the dynamic soil-pile interaction in cohesionless soils by numerical analyses to better characterize the p-y curves considering the nonlinear soil behavior under dynamic loading. A numerical pile-soil-structure model was created in FLAC(3D) and verified by two centrifuge tests published in the literature. The parametric analyses were performed to obtain the p-y curves for various pile diameters, soil relative densities, and degrees of nonlinearities. Based on the parametric analyses, a mathematical model was proposed for the dynamic p-y curves for cohesionless soils. The proposed model characterizes the backbone of dynamic p-y curves based on the three leading parameters (initial stiffness K-py, ultimate resistance p(u), and degree of nonlinearity n). The numerical analyses showed that the p-y curve nonlinearity mainly depends on the employed modulus reduction curves of soils. In the model, the degree of nonlinearity parameter (n) was directly related to the soil parameter gamma(r) which solely represents the modulus reduction curve of soils. In this regard, by correlating the dynamic p-y curves to the reference strain, the dependence on various dynamic soil parameters was diminished. The validation analyses performed in structural analysis software demonstrated that the proposed dynamic p-y model is capable of estimating the pile and structure response under earthquake loading more accurately by incorporating the hysteretic nonlinear soil behavior.