This paper presents a modelling technique for the nonlinear analysis of composite steel-concrete beams with partial shear interaction. It extends the applicability of two stiffness elements previously derived by the authors using the direct stiffness method, i.e. the 6DOF and the 8DOF elements, to account for material nonlinearities. The freedoms are the vertical displacement, the rotation and the slip at both ends for the 6DOF stiffness element, as well as the axial displacement at the level of the reference axis for the 8DOF stiffness element. The solution iterative scheme is based on the secant method, with the convergence criteria relying on the ratios of the Euclidean norms of both forces and displacements. The advantage of the approach is that the displacement and force fields of the stiffness elements are extremely rich as they correspond to those required by the analytical solution of the elastic partial interaction problem, thereby producing a robust numerical technique. Experimental results available in the literature are used to validate the finite element proposed in the paper. For this purpose, those reported by Chapman and Balakrishnan (1964), Fabbrocino et al. (1998, 1999) and Ansourian (1981) are utilised; these consist of six simply supported beams with a point load applied at mid-span inducing positive bending moment in the beams, three simply supported beams with a point load applied at mid-span inducing negative bending moment in the beams, and six two-span continuous composite beams respectively. Based on these comparisons, a preferred degree of discretisation suitable for the proposed modelling technique expressed as a function of the ratio between the element length and depth is proposed, as is the number of Gauss stations needed. This allows for accurate prediction of the nonlinear response of composite beams.