What does Born Green equation signify physically? What does Born Green equation obtained from YBG hierarchy for the equilibrium particle densities signify? I mean how can you model the equation into a physical problem?I understood the steps involved in the derivation of the expression, but am still unsure if I understood it physically, am unable to explain it properly.
The book is this.
 A: I'm not sure if there is much to physically understand in the equation itself, the derivation is where all the physical insight takes place (Kirkwood's superposition approximation). 
What of the equation, then? Why is it important and why should anyone care? It's because the BBGKY hierarchy, while exact, cannot be solved: It is a relation between the (reduced phase-space) distribution functions $f^{(n)} = \mathcal{F} (f^{(n+1)})$, and somehow this infinite dependence must be brought to an end if we wish to solve for anything. Making assumptions about the matter, Born and Green wrote a closure, $f^{(3)} = \mathcal{G}(f^{(2)})$, and now that $f^{(2)} = \mathcal{F} (f^{(3)}) = \mathcal{F} (\mathcal{G}(f^{(2)}))$, we have an equation for $f^{(2)}$ that we can in principle solve. Using the BBGKY hierarchy on $f^{(2)}$, then, we get all the other distribution functions. 
So, moving on to the subject matter, the Born-Green equation reads
$$-k_BT\nabla_1(\log g(r_{12}) + \beta v(r_{12})) = \rho\int\nabla_1v(r_{13})g(r_{13})(g(r_{23})-1)\mathrm{d}\mathbf{r}_3$$
where I've made the simplifying assumption that the potential $v$ and the radial density distribution function $g$ only depend on the distance of the two particles. So given the potential $v$, we can solve this for the structure of the fluid, i.e. for $g$ (from which all the thermodynamic properties can then be calculated). The latter can in fact be obtained very simply from the structure factor, which in turn is the result of SAXS and similar experimental techniques ($g$ can be directly, and trivially, deduced from molecular simulations).
Hill in Statistical Mechanics goes through the motions of solving this analytically using a couple of simplifying approximations. The basic idea is to do the expansion in density $g(r) = e^{-\beta v(r)}(1+\rho g_1(r))$, where $g_1(r)$ is some function (which can be related to the third virial coefficient). Substituting this along with the hard sphere potential (say) to the Born-Green equation, we finally get:
$$g(r) = \left\{\begin{array}{cc}0 , & r \leq a \\ 1 + \frac{4\pi}{3}(\rho a^3)\left(1-\frac{3}{4}\left(\frac{r}{a}\right)+\frac{1}{16}\left(\frac{r}{a}\right)^3\right), & a < r < 2a \\ 
1, & a \leq r\end{array}\right.$$

Finally it should be mentioned that Born-Green is not a particularly good approach: the superposition approximation breaks down quite fast (the book by Hill I mentioned earlier goes on about this for several pages). Other, similar, methods work better. 
Addendum
I decided to solve the equation without linearizing it (i.e. without the dilute approximation). This requires a numerical method, so I tried what first came to mind, fixed-point iteration, and it seems to have worked. First, I recast the integral into the form (or more accurately, I lifted this form from Hill)
$$k_BT \log g(r) = -v(r) + \pi\rho \int_0^\infty v'(s)g(s)\int_ {-s}^s (s^2-y^2) \frac{y+r}{r} (g(y+r) - 1) \mathrm{d}y\mathrm{d}s$$
I massaged this a tad by throwing in the hard sphere potential. I engauged the data off SklogWiki for densities 0.2 and 0.6 and here's the results:

The blue line is 0.6 from SklogWiki, the green the B-G prediction. Red is data for 0.2, and cyan the B-G prediction. 
Finally, the (horrible-looking-but-hopefully-understandable) code I wrote to produce the graphs (I'm not promising it's correct):
L  = 10.
rho = .2
g  = concatenate((zeros(128), ones(128*9)))
dx = 1.*L/(len(g)-1)

r = dx*arange(len(g))
y = arange(-128, 129)*dx
gn = g.copy()

for j in range(300):
  g = .95*g + .05*gn
  for i in range(len(g))[128:]:
    if (i >= len(g)-128):
      kk = i - (len(g)-128) + 1
      grng = concatenate((g[arange(-128, 129-kk) + i], ones(kk)))
    else:
      grng = g[arange(-128, 129) + i]
    gn[i] = exp( -rho*pi*g[128]*sum((1.-y**2) * (y+r[i])/r[i] * (grng-1) * dx) )

I used Python with the roots of SciPy and NumPy imported into the main namespace (or whatever it is that ipython2 --pylab 'qt4' does on my computer).
