A multiscale numerical approach is established to model damage in random glass fiber composites. A representative volume element of a random glass fiber composite is employed to analyze microscale damage mechanisms, such as matrix cracking and fiber-matrix interfacial debonding, while the associated damage variables are defined and applied in a mesoscale stiffness reduction law. The macroscopic response of the homogenized mesoscale damage model is investigated using finite element analysis and validated through experiments. A case study of a random glass fiber composite plate containing a central hole subjected to tensile loading is performed to illustrate the applicability of the multiscale damage model.