Is there a statistical mechanical model with the same qualitative behaviour as polarity? Every chef and chemist knows that oil and water do not mix and that this is due to the polar nature of the water molecule and the non-polar nature of the oil molecule. In fact if one starts the oil in small drops they assemble into larger drops. Is there a natural statistical mechanical model which  has similar predictions? I.e. some number preserving evolution of cells on a grid that exhibit clustering?
 A: One model in which such a behavior can be described precisely is the Ising model.
This can be done in any dimension, but the results are easier to state (and more precise) in the two-dimensional case, so I'll stick to that.
Consider an Ising model in a finite domain in $\mathbb{Z}^2$, say $V_L = \{1,\dots,L\}^2$ with $+$ boundary condition, below the critical temperature $T_c$. Denotes by $m^*(T)>0$ be the spontaneous magnetization.
Let us fix $m\in (-m^*(T),m^*(T))$ and consider the conditional probability measure
$$
\mu^+_{T,L} \Bigl(\cdot \Bigm| \sum_{i\in V_L} \sigma_i = m |V_L|\Bigr),
$$
where $\mu^+_{T,L}$ is the Gibbs measure for the Ising model in $V_L$ at temperature $T$ and with $+$ boundary condition and $|V_L| = L^2$ is the area of $V_L$. I have assumed that $m$ has been chosen such that $m|V_L|$ is a possible value for the magnetization in $V_L$. Note that this conditional probability measure corresponds precisely to the canonical ensemble in the lattice gas language (that is, if we identify $-$ spins with particles and $+$ spins with vacancies, then the measure above is the canonical Gibbs measure with total number of particles given by $\tfrac12(1-m)|V_L|$).
One can then prove that, with probability going to $1$ as $L\to\infty$, configurations sampled from the above measure will have the following structure:

*

*There is exactly one droplet of $-$ phase of diameter of order $L$, immersed inside the $+$ phase;

*Once rescaled by $1/L$, the shape of this droplet converges to a deterministic shape, given by the solution of the following variational problem: minimize the surface tension among all subsets of $\mathbb{R}^2$ with area $\frac{m^*-m}{2m^*(T)}$. The minimizing shape is known as the Wulff shape.

*All other droplets are microscopic (of diameter at most $\log L$).

Since a picture is worth a thousand words, here is a typical configuration with $L=998$ and $m=0$; black pixels correspond to $-$ spins (particles), white pixels to $+$ spins (vacancies):

The restriction to $+$ boundary condition is not a requirement. But considering more general boundary conditions (or even boundary fields) leads to additional interesting phenomena (for instance, the droplet might attach to the boundary, a manifestation of the wetting transition in the Ising model).
If you want additional information on the above, I suggest that you read Section 4.12.1 in this book, or this review paper (more technical, but much more detailed).
