In this thesis, the Monte Carlo simulation method (MCSM) is used in conjunction with finite elements (FEs) for probabilistic analysis of both self-excited vibration and dynamic stability of pipes conveying fluid. For fluid-structure interaction, the Euler-Bernoulli beam kinematic model is used for analyzing pipe structure and plug flow model for representing internal fluid flow through the pipe. Moreover, due to the inherent uncertainty in the system parameters and for determining the problem more precisely, parameters such as masses per unit length of structure and fluid, bending stiffness of structure, and fluid flow velocity are modeled as random fields. By considering the structural and fluid parameters of system as random fields, the deterministic governing partial differential equation (PDE) of continuous system is transformed into a stochastic PDE. The continuous random fields are discretized by mid-point and local average discretization methods. In other words, at the beginning, covariance matrix is generated for each of the random fields by discretization methods. Then, each of those random fields is converted into a random vector by generating an independent Gaussian random vector and using the Cholesky decomposition of covariance matrix. On the other hand, a certain value is assigned to each of the random fields in every finite element due to assuming that the number of FEs is the same as the number of elements required for the discretized random fields. Then, by using the Monte Carlo simulations in each iteration loop, every distributed-parameter PDE having stochastic lumped-parameters is transformed into a deterministic distributed-parameter PDE. Each PDE is transformed into a system of deterministic ordinary differential equations (ODEs) by using FEs. Accordingly, all of the deterministic and stochastic parameters of the system are discretized. For self-excited vibration analysis, the eigenvalue problem is solved for investigating the complex-valued eigenvalues and critical eigenfrequencies. Consequently, having complex eigenfrequencies and divergence velocities, the statistical responses of stochastic problem are obtained like expected values, standard deviations, probability density functions, and the probability of occurrence for divergence instabilities. As expected, the standard deviation values for mid-point discretization method are larger than those for local average. As a result, the mid-point discretization method seems to show the upper bound, while the local average demonstrates to be a lower bound for the standard deviations of the system responses. Furthermore, the probability density function of imaginary part of the eigenfrequency in the first mode appears to consist of two distinct parts: the first part is similar to a Dirac delta function and the second part is similar to a continuous density function. Actually, the former represents the unstable region, while the latter represents the stable part of the system. On the other hand, it has been shown that, the randomness effects of the fluid parameters on the system are much more pronounced than those of structural parameters. Consequently, randomness in the structural parameters could be ignored compared to that in the fluid parameters Keywords: Random vibration, Random fields, Monte Carlo simulation, Fluid-structure interaction, Stochastic stability analysis, Pipes conveying fluid