Satellite networks are promising to provide ubiquitous and high-capacity global wireless connectivity. Traditionally, satellite networks are modeled by placing satellites on a grid of multiple circular orbit geometries. Such a network model, however, requires intricate system-level simulations to evaluate coverage performance, and analytical understanding of the satellite network is limited. Continuing the success of stochastic geometry in a tractable analysis for terrestrial networks, in this paper, we develop novel models that are tractable for the coverage analysis of satellite networks using stochastic geometry. By modeling the locations of satellites and users using Poisson point processes on the surfaces of concentric spheres, we characterize analytical expressions for the coverage probability of a typical downlink user as a function of relevant parameters, including path-loss exponent, satellite height, density, and Nakagami fading parameter. Then, we also derive a tight lower bound of the coverage probability in tractable expression while keeping full generality. Leveraging the derived expression, we identify the optimal density of satellites in terms of the height and the path-loss exponent. Our key finding is that the optimal average number of satellites decreases logarithmically with the satellite height to maximize the coverage performance. Simulation results verify the exactness of the derived expressions.