We apply imaginary-time evolution with the operator to study relaxation dynamics of gapless quantum antiferromagnets described by the spin-rotation-invariant Heisenberg Hamiltonian . Using quantum Monte Carlo simulations to obtain unbiased results, we propagate an initial state with maximal order parameter (the staggered magnetization) in the spin direction and monitor the expectation value as a function of imaginary time . Results for different system sizes (lengths) exhibit an initial essentially size independent relaxation of toward its value in the infinite-size spontaneously symmetry broken state, followed by a strongly size dependent final decay to zero when the rotational symmetry of the order parameter is restored. We develop a generic finite-size scaling theory that shows the relaxation time diverges asymptotically as , where is the dynamic exponent of the low-energy excitations. We use the scaling theory to develop a practical way of extracting the dynamic exponent from the numerical finite-size data, systematically eliminating scaling corrections. We apply the method to spin- Heisenberg antiferromagnets on two different lattice geometries: the standard two-dimensional (2D) square lattice and a site-diluted 2D square lattice at the percolation threshold. In the 2D case we obtain , which is consistent with the known value , while for the site-diluted lattice we find or , where is the fractal dimensionality of the percolating system. This is an improvement on previous estimates of . The scaling results also show a fundamental difference between the two cases; for the 2D square lattice, the data can be collapsed onto a common scaling function even when is relatively large, reflecting the Anderson tower of quantum rotor states with a common dynamic exponent . For the diluted 2D square lattice, the scaling works well only for small , indicating a mixture of different relaxation-time scalings between the low-energy states. Nevertheless, the low-energy dynamic here also corresponds to a tower of excitations.