The Eulerian k-fluid CFD model was used to simulate the macroscale multiphase flow in packed beds. The geometric complexity of the bed structure is resolved by statistically describing the porosity distribution. The complicated multiphase interactions are computed using the Ergun type of formula developed based on bench-scale hydrodynamic experiments. The work is presented in two sequential articles. Part I discusses implementation issues of the k-fluid CFD model for packed beds. The drag exchange coefficients are obtained from the model of Holub et al. for the particle-fluid interfaces X-ks and from the model of Attou et al. (1999) for the gas-liquid interface, X-gl. The effect of particle external wetting on flow distribution was incorporated into the model through the capillary pressure evaluated by either the J-function of Leverett (1941) for air-water or by the expression of Attou and Ferschneider (1999) for other fluids. In the framework of CFDLIB, the choice of the grid size and boundary conditions are discussed. An appropriate relationship between the section size and valiance of the sectional porosity distribution was used for flow simulation. Part 11 discusses the extensive numerical results, and the CFD model is compared with experimental data in the literature.