The Preconditioned Conjugate Gradient (PCG) method is one of the most popular iterative methods for solving large linear systems with a symmetric and positive semi-definite coefficient matrix. However, if the preconditioned coefficient matrix is ill-conditioned, the convergence of the PCG method typically deteriorates. Instead, a two-level PCG method can be used. The corresponding two-level preconditioner usually treats unfavorable eigenvalues of the coefficient matrix effectively, so that the two-level PCG method is expected to converge faster than the original PCG method. Many two-level preconditioners are known in the fields of deflation, multigrid and domain decomposition methods. Several of them are discussed in this thesis, where the main focus is on the deflation method. We show some theoretical properties of the deflation method, which give insights into the effectiveness of this method. A crucial component of the deflation preconditioner is the choice of projection vectors. Several choices are discussed and examined. We advocate that subdomain projection vectors, which are based on disjoint and piecewise-constant vectors, are among the best choices for a class of problems. Subsequently, we examine the application of the deflation method to linear systems with singular coefficient matrices. Several mathematically equivalent variants of the original deflation method are proposed to deal with the possible singularity of this coefficient matrix. In addition, two approaches are discussed in order to handle coarse linear systems with a Galerkin matrix, which are involved in each iteration of the deflation method. After the discussion of the implementation and efficiency issues of the deflation method, it is demonstrated that this method is usually faster than the original PCG method. Moreover, we present a comparison between the deflation method and other well-known two-level PCG methods, among them the balancing-Neumann-Neumann, additive coarse-grid correction, and multigrid methods based on symmetric and nonsymmetric V-cycles. As the parameters of the corresponding two-level preconditioners are abstract, we show that these methods are strongly connected to each other. The comparison is also done where the different two-level PCG methods adopt their typical and optimized set of parameters. Numerical experiments show that some multigrid methods are attractive in addition to the deflation method. The major application of this thesis is the Poisson equation with a discontinuous coefficient, which is derived from 2-D and 3-D bubbly flow problems. Most of the performed numerical experiments in this thesis are based on this equation. Both stationary and time-dependent experiments are carried out to emphasize the theoretical results. We show that two-level PCG methods are significantly faster than the original PCG method in almost all experiments. Hence, computations involved in bubbly flows can be performed very efficiently using these PCG methods.