  In this paper, we use large $\pppm$ N-body simulations to study the
three-point correlation function $\zeta$ of clusters in two theoretical models.
The first model (LCDM) is a low-density flat model of $\Omega_0=0.3$,
$\Lambda_0=0.7$ and $h=0.75$, and the second model (PIM) is an
Einstein-de-Sitter model with its linear power spectrum obtained from
observations. We found that the scaled function $Q(r,u,v)$, which is defined as
the ratio of $\zeta (r, ru, ru+rv)$ to the hierarchical sum $\xi (r)\xi (ru)+
\xi (ru) \xi (ru+rv) +\xi (ru+rv)\xi (r)$ (where $\xi$ is the two-point
correlation function of clusters), depends weakly on $r$ and $u$, but very
strongly on $v$. $Q(r,u,v)$ is about 0.2 at $v=0.1$ and 1.8 at $v=0.9$. A model
of $Q(r,u,v)=\Theta 10^{1.3v^2}$ can fit the data of $\zeta$ very nicely with
$\Theta\approx 0.14$. This model is found to be universal for the LCDM clusters
and for the PIM clusters. Furthermore, $Q(r,u,v)$ is found to be insensitive to
the cluster richness. We have compared our N-body results with simple
analytical theories of cluster formation, like the peak theories or the local
maxima theories. We found that these theories do not provide an adequate
description for the three-point function of clusters. We have also examined the
observational data of $\zeta$ presently available, and have not found any
contradiction between the observations and our model predictions. The
$v$-dependence of $q$ in a projected catalogue of clusters is shown to be much
weaker than the $v$-dependence of $Q$ found in the three-dimensional case. It
would be important to search for the $v$-dependence of $Q$ in redshift samples
of rich clusters.
