The Discrete Element Method (DEM) has been used to model elastic, long fibers in diverse areas [1]. Initially, the fibers were discretized with spheres, but the need for a higher level of accuracy and less computing cost led to the use of non-spherical and elongated particles. The contact model presented here employs sphero-cylinders and was specifically developed to compute the dynamics of curved fibers with heterogenous radii. It differentiates between shear and axial forces as well as bending and torsion torques and can predict possible fiber failure under different scenarios. The disruptive idea of this model is the presence of two springs (instead of one) for each force and torque considered in the bond, one at each side of the bond, without the need to average the normal or tangential directions. The stiffness of each spring can be computed from the physical properties, the cross-section, and the length of the fibres of each side, as partially explained in [2]. The decomposition of the angle into bending and torsion must be done in a univocal way, avoiding popular non-univocal decompositions like the Euler angles. This model uses the bi-normal vector to decompose the total relative angle into one bending angle around the bi-normal vector and one torsion angle. Because of the possible heterogeneity of the particles and the different nature of shear and axial stiffness (or bending and torsional stiffness), this model accepts that the stiffness of the springs at both sides of the bond is different, so the bond itself needs to be re-located and re-oriented at every time step to ensure that the forces and torques generated by the springs at both sides are equal. The limits of the elastic range can be defined in different ways (stress, angle). Once reached, depending on the simulated material, other failure modes can be used: plasticity, breakage (bond rupture), or buckling (reversible buckling, damage buckling, plastic buckling).