A three-dimensional, implicit, finite element composites damage model is developed, and applied to a problem involving highly complex, three-dimensional loading, i.e. a single-lap, multi-bolt, composite joint, having variable clearances. The model extends a previous model, adding Puck's criteria to allow delamination-like, transverse cracking-like or mixed mode failures to be simulated. Plasticity is also added for non-reversible strains in the undamaged parts of the matrix, and a mesh independence strategy is incorporated. A back-up strategy for when the Newton Raphson method fails to converge, adds significantly to model robustness. The model is implemented in a commercial, implicit solver and demonstrates excellent robustness, being capable of following damage progression from onset, through bearing failure at the holes, to catastrophic, net-tension failure. The predicted failure modes and loads are in good agreement with experiment The effect of clearance on secondary bending, and consequently on damage progression, is demonstrated. The current model shows better capability than our previous model for predicting matrix damage near the hole, which will be important for future fatigue modelling studies. (C) 2015 Elsevier Ltd. All rights reserved.