What is the statistical equilibrium for this simulation of happy bubbles? Happy Bubbles

I hope this is not too specific or practical, but I just made a simulation of sorts and seem to have hit quite close to an equilibrium (by accident). Now I am wondering if and how you would calculate the exact function or exact value for this situations statistical "sweetspot".

The Problem

I was playing around with D3.js and created a canvas with randomly appearing bubbles with a 5% chance every iteration. Also every iteration bubbles touching any border will disappear. This can happen because the bubbles move around every iteration.

The Parameters


Starting with a single bubble
Every bubble has a x- and a y-coordinate which represents the center of the bubble
Bubble Size: 20px radius
Iterationlength: 100 ms
Canvas Size: 600px by 600px
Bubble Speed: -5px to 5px on X-Axis and Y-Axis

(Math.random() * 10) - 5 (the random number is a floating point number between 0 (including) and 1 (excluding))

Starting position per bubble: Math.random() * 600 for x and y coordinate (seperately)
Bubble pops when it gets within 20px of any outer edge (again using the center as reference)


The Equilibrium

With these Parameters new bubbles and disappearing bubbles seem to cancel each other out. The number of bubbles ranges from 100 to 130 but has not breached those limits in more then an hour.

The Question

Do these Parameters ensure such an equilibirum or is it nothing special, that the number seems to be rather stable? How would one calculate this?

Caveats

One I can think of is that the bubble could move from 20.000001 to 15.000001 in one iteration. But I guess that is a small detail. Maybe this isn't even relevant.

Edit:

Here is the jsFiddle! (Be carefull with less then 100 ms iteration time, it might stress the browser)
This is a simple model to estimate $\mathcal{N}$, the expected number of bubbles after some sort of equilibrium is reached. In this model, the bubble will be randomly displaced at each iteration.

The first claim is $\mathcal{N}$ is proportional to $\mathcal{T}$, the expected life time of a bubble . After we start the system and wait long enough for the system to stabilize, there are two ways to estimate $\mathcal{N}$.


Observe the system for $N \gg 1 $ iterations. Sum the number of bubbles observed in these iterations and then divide the final result by $N$.
Follow all bubbles created and destroyed during these $N$ iterations. Sum over their life time and then divide by $N$.


In general, these two ways of counting will differ by a number which is of the order of $\mathcal{T}$. If $\mathcal{T}$ is finite.
This two method will produce the same estimation of $\mathcal{N}$ as $N \to \infty$.

If on average, $\alpha$ bubble is created during each iteration, one will have

$$\mathcal{N} = \alpha \mathcal{T}$$

In this problem $\alpha = 0.05$. If the bubble is created within $20$px from the edge, the bubble will be destroyed in the next iteration. 
The same thing happens to the bubble at the end of its life. All this process
will add $1$ to the expected life time $\mathcal{T}$.

For a probability $\left(\frac{560}{600}\right)^2$, the bubble will be created
within the safe zone $[20,580] \times [20,580]$.
Let $L = 580-20 = 560$ and $t$ be the number of iteration since the bubble 
is created. Introduce random variables $X, Y, \eta_x(t), \eta_y(t)$ and we will
update the position of the bubble by following recipe

$$\begin{cases}
x(0) &amp;= 20+X\\
y(0) &amp;= 20+Y
\end{cases}
\quad\text{ and }\quad
\begin{cases}
x(t) &amp;= x(t-1) + \eta_x(t)\\
y(t) &amp;= y(t-1) + \eta_y(t)\\
\end{cases}
\quad\text{ for }\quad t &gt; 0
\tag{*1}
$$

Subject to the constraint the bubble is created in the safe-zone, the variables $X, Y$ and $\eta_x(t), \eta_y(t)$ are uniform random variables taking values over $[0,L]$ and $[-5,5]$ respectively.

Let $T$ be the random variable corresponds to the life time of a bubble (excluding the last iteration discussed above).
Introduce a function
$$u(x,y) = {\bf E}\big[ T \big| (X,Y) = (x,y) \big]\quad\text{ for }\quad (x,y) \in [0,L]\times[0,L]$$
If $u(x,y)$ can be approximated by a continuous function, then the 
$\mathcal{T}$ we want will be given by a integral

$$\mathcal{T} = 1 + \left(\frac{560}{600}\right)^2 {\bf E}[T]
\quad\text{ with }\quad
{\bf E}[T] = 
\frac{1}{L^2} \int_0^L\int_0^L u(x,y) dx dy
$$

We can view the random update in $(*1)$ as a Brownian motion. If equilibrium has been reached, we should have

$$u(x,y) = {\bf E}\big[ u(x + \eta_x, y + \eta_y) \big] + 1$$
This implies $u(x,y)$ should approximately satisfies a PDE

$$u(x,y) = u(x,y) 
+  \frac{{\bf E}[\eta_x^2]}{2} \frac{\partial^2 u(x,y)}{\partial x^2}
+  \frac{{\bf E}[\eta_y^2]}{2} \frac{\partial^2 u(x,y)}{\partial y^2}
+ 1
$$
Let $\sigma^2 = {\bf E}[\eta_x^2] = \frac{5^2}{3}$, this has the form of a Possion
equation
$$\nabla^2 u(x,y) = -\frac{2}{\sigma^2}$$
and the appropriate condition for this problem is $u(x,y) = 0$ on the boundaries. 

If I'm not mistaken, the solution of these boundary value problem is given by

$$u(x,y) = \frac{32L^2}{\sigma^2\pi^4} \sum_{n=1,\text{odd}}^\infty\sum_{m=1,\text{odd}}^{\infty} 
\frac{1}{n m (n^2+m^2)} 
\sin\left(\frac{\pi n x}{L}\right)
\sin\left(\frac{\pi m y}{L}\right)
$$
Integrate this, we get

$$\mathcal{T} = 1 + \left(\frac{560}{600}\right)^2\frac{128L^2}{\sigma^2\pi^6} K
\quad\text{ where }\quad
K = \sum_{n=1,\text{odd}}^\infty\sum_{m=1,\text{odd}}^{\infty} 
\frac{1}{n^2 m^2 (n^2+m^2)} 
\approx 0.52973
$$
Substitute this into the formula $\mathcal{N} = \alpha T$, we get $$\mathcal{N}_{600} \approx 115.65$$
which now falls within the observed range of OP's numerical result. 
If one repeat the same argument to other canvas size $s$, we have

$$\mathcal{N}_s = 1 + \left(1 - \frac{40}{s}\right)^2
\frac{384(s - 40)^2}{25\pi^6} K
$$
Following is a table for $\mathcal{N}_s$ for other canvas sizes $s$:

$$
\begin{array}{r:r}
s &amp; \mathcal{N}_s\\
\hline
300 &amp;  21.54\\
400 &amp;  44.47\\
500 &amp;  75.84\\
600 &amp; 115.65\\
700 &amp; 163.92\\
800 &amp; 220.64\\
900 &amp; 285.82
\end{array}
$$
I fitted the average path lengths I obtained.
$x$ axis: Size of square canvas $\times 100$ pixels.
$y$ axis: Average number of iterations per bubble.