A numerical method was developed to simulate the fouling processes on tubes. The detailed particle deposition and removal mechanisms were included in the model and the evolution of the shape of fouling layers was obtained. Multiple-relaxation-time lattice Boltzmann method (MRT-LBM) and finite volume method (FVM) were coupled to simulate the air flow. The particle motion was simulated by the probabilistic cellular automata model. The restitution coefficient was calculated by energy conservation and used as the deposition criterion. The particle removal was determined by the force and moment analysis. Since the simulation time was much shorter than the real time of the fouling process, a ratio was proposed for the time conversion between the simulation time and real time. The fouling processes for different particle diameters and inlet velocities were simulated by the proposed method. When the mass concentration was specified, small particles had large fouling rates. The fouling area grew linearly with time without removal mechanism, but grew exponentially to an asymptotic balance value when the removal mechanism was considered. The mass flux, particle deposition and removal were simultaneously influenced by the inlet velocity, so the relation between the inlet velocity and fouling rate was not monotonic. A velocity range existed in which the fouling rate was high. The removal was severe on the windward side and the area right behind the tubes, but relatively moderate on the laterals of the leeward side. The fouling layers grew on the entire leeward side of the tubes. On the windward side, the cone-shaped fouling layers were formed, which changed the flow and stopped the further particle deposition. (C) 2016 Elsevier Ltd. All rights reserved.