Structural topology optimization with random material or boundary load requires significant computational costs in solving numerous random elastic systems, particularly in 3D. Based on the multimode representation of displacement and the Monte Carlo method for sampling, a phase field model is proposed for robust topology optimization. The well-posedness and truncation error estimates of the multimode Monte Carlo method have been discussed. Numerical examples are presented.