A computational scheme for three-dimensional (3-D), discrete simulation of bubble dynamics in gas-liquid bubble columns and gas-liquid-solid fluidized beds is presented. This scheme describes the motion of the gas, liquid, and solid phases, respectively, based on the level-set interface tracking method, the locally averaged time-dependent Navier-Stokes equations coupled with the Smagorinsky Sub-Grid Scale stress model, and the Lagrangian particle motion equations. Various gas, liquid, and solid interactive forces are considered, which include those due to drag, added mass, intefacial tension, bubble-particle collision and particle-particle collision. The accuracy and validity of the simulation based on this scheme is examined through computational experiments on some well-known bubbling phenomena including the rise behavior of a single gas bubble in liquids and the bubble formation process from an orifice embedded in liquids. The computed bubble rise velocity, bubble shapes and their fluctuations, and bubble formation behavior in liquids show a good agreement with existing theories and experimental findings. The simulation is extended to investigate the effects of solids concentration on bubble formation and rise behavior. It is found that the bubbles are less coalesced in a fluidized bed with a higher solids concentration. A circulatory motion of particles with an upward flow in the center region and a downward flow near the wall, and a lateral solids concentration distribution in a gas-liquid-solid fluidized bed can also be predicted. The bubble formation and rise behavior in a multi-nozzle bubble column are simulated, and the limitation of the Smagorinsky Sub-Grid Scale stress model is discussed (C) 2004 American Institute of Chemical Engineers.