The objective of this work is to develop a microstructure-based simulation approach to assess the fatigue life of solder joints that are used by the microelectronics industry. The developed approach can generate solder joints with random grain morphologies by means of 3D Voronoi tessellation. The anisotropic material behavior of each grain is described by the Garofalo creep equation combined with Hill's definition of the equivalent stress for anisotropic materials. Grain boundaries are implemented as interface elements, with an isotropic creep constitutive model. The stochastic variability in the creep response of solder joints was qualitatively estimated by generating 100 unique solder joints containing 5 to 9 grains, each having a random material orientation. These joints were independently loaded with a realistic stress level for microelectronic products during thermal cycling. The volume-averaged creep strain energy density in the solder joints was used to predict the fatigue life of the solder joints. The results showed a factor of ~4 difference in expected lifetime of the individual solder joints. Next, nine randomly picked solder joints from the above-mentioned pool of 100 were sandwiched between a silicon die and a printed circuit board to form a simulation model of a Wafer-Level Chip-Scale package (WLCSP). The creep strain energy density in the joints was computed for 34 unique cases of the WLCSP. A factor of ~2.5 between the highest and lowest estimate for the solder joint life was found. The slope of the corresponding Weibull distribution equals ~6, which falls within the slopes typical reported for solder joint reliability of WLCSPs.