The paper developed a new numerical joint model based on nonlinear finite element method for the multiple bolted joints in composite laminates. The new joint model consists of the geometrically exact beam and shell elements, as well as the new Hertz contact elements. The geometrically exact elements, beams and shells, are used to model the bolts and laminates with nonlinear deformations, respectively. Meanwhile, the new contact elements apply the Hertz model to simulate the contact, slippage and bolt-laminate interaction. This contact element combined with the multi-point constraints to transfer the shear forces from one shell through beam to another shell element. As an essential part of the bolted joint model, a new type of geometrically exact shell element was developed with second-order accuracy and nine nodes. The strain-displacement relationship of this shell element was formulated according to the rotation-free Green-Lagrange strain tensor. This shell element approximated the higher-order components of strain tensor, ensuring it can accurately describes the real deformation of thin-walled composite laminates. Application of this joint model to a three-bolt, single-shear joint in composite laminates was presented, and its predictions are compared with those of commercial code ABAQUS. Numerical simulation results show that the new joint model is suitable for the design of bolted joints in composite laminates with good accuracy and high efficiency.