A design of experiment (DOE) methodology based on numerical simulation is presented to improve thermal fatigue reliability of multirow quad flat nonlead (QFN) packages. In this method, the influences of material properties, structural geometries, and temperature cycling profiles on thermal fatigue reliability are evaluated, a L27(38) orthogonal array is built based on Taguchi method to figure out optimized factor combination design for promoting thermal fatigue reliability. Analysis of variance (ANOVA) is carried out to examine the influence of factors on the thermal fatigue reliability and to find the significant factors. Anand constitutive model is adopted to describe the viscoplastic behavior of lead-free solder Sn3.0Ag0.5Cu. The stress and strain in solder joints under temperature cycling are studied by 3D finite element (FE) model. The modified Coffin–Manson model is employed to predict the fatigue life of solder joints. Results indicate that the coefficients of thermal expansion (CTE) of printed circuit board (PCB), the height of solder joint, and CTE of epoxy molding compound (EMC) have critical influence on thermal fatigue life of solder joints. The fatigue life of multirow QFN package with original design is 767 cycles, which can be substantially improved by 5.43 times to 4165 cycles after the optimized factor combination design based on the presented method.