A novel strain-displacement variational formulation for the flexural behaviour of laminated composite beams is presented, which accurately predicts three-dimensional stresses, yet is computationally more efficient than 3D finite element models. A global third-order and layer-wise zigzag profile is assumed for the axial deformation field through the laminate thickness to account for the effect of both stress-channelling and stress localisation. By using this axiomatic kinematic field, a variational formulation is developed based on an equivalent single layer theory such that the number of unknowns is independent of the number of composite layers. The axial and couple stresses are evaluated from the displacement field, while the transverse shear and transverse normal stresses are computed by the interlaminarcontinuous equilibrium conditions within the framework of the modified couple stress theory. Axial and transverse force equilibrium conditions are imposed via two Lagrange multipliers, which correspond to the axial and transverse displacements. Using this mixed variational approach, both displacements and strains are treated as unknown quantities, resulting in more functional freedom to minimise the total strain energy. In this work, two strain-displacementbased models are developed to investigate the effect of couple stress on the flexural behaviour of composite beams. The first model neglects the presence of couple stress, whilst the second includes couple stress with an additional unknown curvature. The differential quadrature method is used to solve the resulting governing and boundary equations for simply-supported and clamped laminated beams. For the simply-supported case, numerical results from this variational formulation for both models agree well with those from a Hellinger-Reissner stressdisplacement mixed model found in the literature and the 3D elasticity solution given by Pagano. For the clamped laminate, the additional curvature associated with the couple stress plays an important role in accurately predicting stresses, which is confirmed by a high-fidelity 3D finite element model.