Direct simulations of flexible cylindrical fiber suspensions in a finite Reynolds number flow are reported. The simulation method is based ona lattice Boltzmann equation and a flexible fiber model. A slender solid body is discretized into a chain of cylindrical segments contacting each other at the fiber ends through ball and socket joints that allow adjunct segments rotate around the joints in three dimensional space. A constraint force is imposed at each joint. In general, motion and rotational matrix of each segment are functions of constraint forces. It is necessary to linearize rotational matrix in forces and torques so that constraint forces could be solved using joint contacting conditions. Therefore, quaternion parameters as well as rotational matrix could be expanded in a power series of the length of time step up to a second order. A half frog leap algorithm is modified to ensure the ball and socket joint conditions satisfied at each time step. The validation of the present flexible fiber method is tested by using a rigid particle method when the fiber is very stiff. It is shown that the computational results are consistent with the existed experimental and theoretical results at finite Reynolds number flows. With present method, nonlinear inertial interactions between fluid and flexible filament can be naturally studied.