In this paper, a new finite element model has been introduced with the aim of efficient investigation of residual thermal stresses in fiber-reinforced composites, in which the inhomogeneous interphase is considered. For the inhomogeneous interphase modeling, four different kinds of material properties variation of the interphase (power, reciprocal, cubic and exponential variations) with the radial coordinate have been used. A mono fiber circular unit cell is considered using a finite element (FE) method. Extending the mono fiber model, FE models with different arrays of fibers have been created to investigate the effects of neighboring fibers on the results. In order to assure the convergence of results, a convergence analysis has been carried out for each of the models. To verify the finite element model, the FE results are compared with theoretical results available in the literature. In this paper, three different types of RVE configurations, circular, square and hexagonal are modeled and the effects of each type of fiber packing are studied. Performing an extensive study, the appropriate boundary conditions for RVEs are presented. The boundary conditions presented in this research are proved to be able to model the overall behavior efficiently.