Implicit large eddy simulations of two-dimensional Rayleigh-Taylor instability at different density ratios (i.e., Atwood number A=0.05, 0.5, and 0.9) are conducted to investigate the late-time dynamics of bubbles. To produce a flow field full of bounded, semibounded, and chaotic bubbles, three problems with distinct perturbations are simulated: (I) periodic sinusoidal perturbation, (II) isolated W-shaped perturbation, and (III) random short-wave perturbations. The evolution of height h, velocity v, and diameter D of the (dominant) bubble with time t are formulated and analyzed. In problem I, during the quasisteady stage, the simulations confirm Goncharov's prediction of the terminal speed v_{∞}=Frsqrt[Agλ/(1+A)], where Fr=1/sqrt[3π]. Moreover, the diameter D at this stage is found to be proportional to the initial perturbation wavelength λ as D≈λ. This differed from Daly's simulation result of D=λ(1+A)/2. In problem II, a W-shaped perturbation is designed to produce a bubble environment similar to that of chaotic bubbles in problem III. We obtain a similar terminal speed relationship as above, but Fr is replaced by Fr_{w}≈0.63. In problem III, the simulations show that h grows quadratically with the bubble acceleration constant α≡h/(Agt^{2})≈0.05, and D expands self-similarly with a steady aspect ratio β≡D/h≈(1+A)/2, which differs from existing theories. Therefore, following the mechanism of self-similar growth, we derive a relationship of β=4α(1+A)/Fr_{w}^{2} to relate the evolution of chaotic bubbles in problem III to that of semibounded bubbles in problem II. The validity of this relationship highlights the fact that the dynamics of chaotic bubbles in problem III are similar to the semibounded isolated bubbles in problem II, but not to that of bounded periodic bubbles in problem I.