We explore the signature of phase correlations in Fourier modes of dark matter density fields induced by nonlinear gravitational clustering. We compute the distribution function of the phase sum of the Fourier modes, theta(k1) + theta(k2) + theta(k3), for triangle wavevectors satisfying k(1) + k(2) + k(3) = 0 and compare with the analytic prediction in perturbation theory recently derived by one of the authors of this paper. Using a series of cosmological N-body simulations, we extensively examine the time evolution and the dependence on the configuration of triangles and the sampling volume. Overall, we find that the numerical results are remarkably consistent with the analytic formula from the perturbation theory. Interestingly, the validity of the perturbation theory at a scale k corresponding to the wavevector k is determined by P(k)/V-samp, the ratio of the power spectrum P(k) and the sampling volume V-samp, not by k(3)P(k), as in the case of conventional cosmological perturbation theory. Consequently, these statistics of phase correlations are sensitive to the size of the sampling volume itself. This feature does not show up in more conventional cosmological statistics, including the one-point density distribution function and the two-point correlation functions, except as a sample-to-sample variation. Similarly, if the sampling volume size Vsamp is fixed, the stronger phase correlation emerges first at the wavevector for which P(k) becomes largest, i.e., in linear regimes according to the standard cosmological perturbation theory, while the distribution of the phase sum stays fairly uniform in nonlinear regimes. The above feature can be naturally understood from the corresponding density structures in real space, as we discuss in detail.