We have performed a series of N-body simulations in cold dark matter (CDM) cosmologies in order to study gravitational clustering properties of clusters of galaxies. Specifically, we consider two spatially flat CDM models: (OMEGA, lambda, h, b) = (1.0. 0.0, 0.5, 1.7) and (0.2, 0.8, 0.75, 1.0), where OMEGA, lambda, h, and b are the cosmological density parameter, the dimensionless cosmological constant, the Hubble parameter in units of 100 km s-1 Mpc-', and the bias parameter, respectively. Clusters of galaxies in the simulations are identified applying an adaptive linking method. Their spatial clustering properties are compared with those for initial density peaks which would be possible sites for cluster formation. The latter predictions are independently obtained on the basis of random-Gaussian theory and from direct simulations. Almost all massive clusters identified in the simulations turn out to contain high-density peaks in qualitative agreement with the previously proposed idea of biased structure formation. Two-point correlation functions for clusters are approximately given by a power-law form and significantly exceed underlying mass correlations. The observed richness-dependent enhancement of the cluster correlation is well reproduced in low-density CDM models, while OMEGA = 1 models show weaker cluster correlation compared to the observational counterpart. Three-point correlation function and skewness of the clusters and peaks in these simulations are also discussed.