Steady-state and time-dependent two-dimensional thermal convection in a Boussinesq, infinite-Prandtl-number fluid with stress-free boundaries has been investigated. Two independent numerical methods have been employed to calculate the evolution of convective flows in a rectangular box with aspect ratio λ=1.8 in a Rayleigh-number (Ra) range of 106<Ra<109. With increasing Ra, greater than 107, the flow reveals the presence of disconnected thermals, rather than connected plumes, driven by a persistent large-scale circulation. Such features have also been reported from laboratory convection experiments in the regime of hard turbulence. Extensive calculations were performed (up to 140 overturns) in order to reach the statistically stationary regime for strongly chaotic flows. A Gaussian distribution with a mean value Nut was derived from the time history of the Nusselt (Nu) numbers. The value of Nut can be directly obtained by solving the steady-state equations via an iteration procedure. Thus the stationary flow obtained from the steady-state method resembles the turbulent flow in a statistical sense. Since the iteration procedure is about 104 times faster than calculating the full time-dependent evolution, it allows for the systematic investigation of the heat-transfer Nu-Ra relationship and other types of scaling laws. The steady-state and time-dependent experiments indicate that a power-law exponent of β=0.315 holds for the Nu-Ra relation for stress-free boundaries in the entire range of Ra. No indication of a jump in the exponent was found in the transition to hard turbulence.