This study extends the curvilinear grid finite‐difference method (FDM) to 3D time‐domain seismic wave modeling through the first‐order velocity‐stress wave equations in generally anisotropic media involving a topographic free surface. In this method, collocated grids are utilized to define all the necessary anisotropic wavefield components at the same grid point to avoid additional errors caused by interpolation in the staggered grid FDM and calculation of the Hooke sum. The complex topographic free surface is discretized by curvilinear grids, which properly describe the free surface and avoid the numerical scattering caused by the staircase approximations of the Cartesian grids. The free‐surface boundary condition in generally anisotropic media is implemented by the traction image method, whose formula in this media is also derived. To validate the proposed method, the numerical solutions of this algorithm are compared with a set of totally independent spectral element solutions. The results present a suitable agreement with the reference solutions to confirm the validity of the algorithm. The influences of the irregular free surface on the anisotropic seismic wavefield are also demonstrated in numerical experiments.