Particle assembly process and settling behaviours in it fluid are studied with a small-population particle system involving in interparticle force. Particle motions are solved individually by the equation of linear and angular momentums, and fluid phase is solved by direct numerical simulation resolving from the boundary layer of the particle surface to the far wake. Interaction between the solid and fluid phases is treated by an immersed boundary method proposed by the present authors (Kajishima et al., JSME Int. J. Ser. B 44-4 (2001) 526; Int. J. Heat and Fluid Flow 23 (5) (2002) 639). The particles form agglomerate due to the interparticle force. The descent velocities of the agglomerates of' various sizes and interparticle strengths are higher than that of a single particle. By varying the parameters of the interparticle forces, transitions in shape and descent velocity of the agglomerates are observed depending on the balance of interparticle and hydrodynamic forces. The following correlations are found: fluid drag per unit area of the agglomerate with respect to the area of the convex closure of the agglomerate projected on a plane normal to the descent direction, drag coefficients of the agglomerate with respect to the squared agglomerate Reynolds number defined with descent velocity and equivalent diameter of the convex closure. (C) 2007 Elsevier B.V. All rights reserved.