Using a new, high mass resolution (Δmb = 105.5 M☉) hydrodynamic simulation of a spatially flat ΛCDM cosmological model with detailed microphysics and galaxy formation, including radiation shielding, energy deposition, and metal enrichment from supernovae and associated metal cooling/heating, we compute the metallicity evolution of damped Lyα systems (DLAs) and find a reasonable agreement with observations. In particular, the observed slow evolution of the DLA metallicity occurs naturally in the simulation as a result of the combined effects of physical and observational selection. The slow metallicity evolution is caused by the steady transformation, with increasing time, of the highest metallicity systems to galaxies, thus depleting this category, while all the lower metallicity systems show, individually, an increase in metallicity. Although the trend of DLA metallicity with redshift is in good agreement with observations, it appears that the average metallicity of simulated DLAs is higher than observed by 0.3-0.5 dex in the probed redshift range (z = 0-5). Our study indicates that this difference may be attributed to observational selection effects due to dust obscuration. If we allow for a dust obscuration effect, our model reproduces the observed metallicity evolution in both amplitude and slope. We find that DLAs are not a simple population but probe a range of different systems, and the mix changes with redshift. The median luminosity of a DLA, LDLA(z), in units of typical galaxy luminosity at that redshift, L*(z), that is, z, decreases from 1.1 to 0.5 as redshift declines from z = 3 to 0, but the absolute luminosity of the median DLA system increases in the same interval by a factor of 5 from 0.1L*(z = 0) to 0.5L*(z = 0). About 50% of all metals in the gaseous phase is in DLAs at all times from z = 5 to z = 1, making a rapid downturn at z ≤ 1 to ~20% by z = 0, as metals are swept into the hotter components of the intergalactic medium (IGM) as well as locked up in stars. While not the primary focus of this study, we find that the model provides good matches to observations with respect to column density distribution and evolution of neutral gas content, if the same dust obscuration is taken into account. We find ΩDLA,comp = 1-3 × 10-3, depending on the effect of dust obscuration.