A DOE (Design of Experiment) methodology based on finite element analysis is presented to investigate thermal fatigue reliability of multi-row QFN packages. In this method, the influences of material properties, structural geometries and temperature cycling profiles on thermal fatigue reliability are evaluated, a L27(3 8 ) orthogonal array is built based on Taguchi method to figure out optimized factor combination design for promoting thermal fatigue reliability. 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 model. The modified Coffin-Manson model is employed to predict the fatigue life of solder joints. Results indicate that CTE of PCB board, the height of solder joint and CTE of epoxy molding compound have critical influence on thermal fatigue life of solder joints. The fatigue life of multi-row QFN package with original design is 767 cycles, which can be substantially improved by 5.43 times to 4165 cycles after optimized factor combination design.