The Karhunen–Loève (K–L) expansion method is a powerful tool for simulating stationary and nonstationary, Gaussian and non-Gaussian stochastic processes with explicitly known covariance functions. Since the K–L expansion requires the solution of Fredholm integral equation of the second kind, it is generally not feasible to simulate multi-dimensional random fields. This is because even the numerical solution of the Fredholm integral multi-dimensional eigenvalue problem is difficult to obtain. In order to address this problem, this paper develops a consistent generalization of K–L expansion for multi-dimensional random field simulation. The new method decomposes an n-dimensional random field into a total of n stochastic processes, each can be represented by using the traditional K–L expansion. Thus, the developed method is embedded into the well-established framework of the K–L expansion for simulating stochastic process, and obviating the need for solving the multi-dimensional integral eigenvalue problems. Four examples, including random fields with different kinds of covariance functions, are used to demonstrate the application of the proposed method.