A computational scheme for discrete-phase simulation of a gas-liquid-solid fluidization system and a two-dimensional code based on it are developed in this study. In this scheme, the volume-averaged method, the dispersed particle method, and the volume-of-fluid (VOF) method are used to account for the flow of liquid, solid particles, and gas bubbles respectively. The gas-liquid interfacial mass, momentum and energy transfer is described by a continuum surface force (CSF) model. A close-distance interaction (CDI) model is introduced which illustrates the motion of the particle prior to its collision: upon collision, the hard sphere model is employed. The particle-bubble interaction is formulated by incorporating the surface tension force in the equation of motion of particles. The particle-liquid interaction is brought into the liquid phase Navier-Stokes (N-S) equations through the use of Newton's third law of motion. The volume-averaged liquid phase N-S equations are solved using the time-split two-step projection method. The simulation results using this scheme are verified for bed expansion and pressure drop in liquid-solid fluidized beds. The simulation of a single bubble rising in a liquid-solid suspension and the particle entrainment by a bubble on the surface of the bed is conducted and the results are in agreement with the experimental findings.