Based on Euler-Lagrange frame, a true three-dimensional numerical simulation of bubbling fluidized bed embedded with two immersed tubes is presented. The solid phase is composed of 178,200 particles of 883 pm diameter and simulated by discrete element method (DEM, a soft-sphere approach). The gas phase is computed through solving the volume-averaged four-way coupling Navier-Stokes equations in the Smagorinsky SGS tensor model is used in large eddy simulation (LES). Particle-tube collision is particularly treated as a transformation of DEM. The volume segmentation of a particle sphere for fraction calculation is solved via a numerical sub-division approach. The numerical results are compared with the experimental results for validation. The results obtained with and without the LES model are also compared. The numerical results show a strong correlation between gas-particle interaction, particle-particle interaction, pressure drop, particle back mixing motion and bubble motion, and of them follow a similar pattern of synchronous periodic variation though the periodicity may vary depending on different flow conditions. The effects of SGS tensor on evolution of fluidized bed are found in various aspects. Finally, the distribution of particle-tube impact frequency is given. (c) 2008 Elsevier Ltd. All rights reserved.