This paper focuses on the continuum scale modeling of dislocation-grain boundary interactions and enriches a particular strain gradient crystal plasticity formulation (convex counter-part of Yal double dagger inkaya et al., J Mech Phys Solids 59:1-17, 2011; Int J Solids Struct 49:2625-2636, 2012) by incorporating explicitly the effect of grain boundaries on the plastic slip evolution. Within the framework of continuum thermodynamics, a consistent extension of the model is presented and a potential type non-dissipative grain boundary description in terms of grain boundary Burgers tensor (see e.g. Gurtin, J Mech Phys Solids 56:640-662, 2008) is proposed. A fully coupled finite element solution algorithm is built-up in which both the displacement and plastic slips are considered as primary variables. For the treatment of grain boundaries within the solution algorithm, an interface element is formulated. The proposed formulation is capable of capturing the effect of misorientation of neighboring grains and the orientation of the grain boundaries on slip evolution in a natural way, as demonstrated by bicrystal specimen examples.