In non linear concrete structure modeling, an efficient numerical resolution scheme is as important as thechoice of the constitutive law. When using Newton-Raphson approach, providing tangent stiffness matrix iscrucial to improve the numerical efficiency related to the iterative process. However, non-differentiablerelations are often encountered (e.g. the Macauley bracket or absolute operator) does not allow the analyticalderivation of the tangent matrix.In this study, three strategies are proposed and compared to deal with this issue. These strategies are appliedto the Mazars Unilateral damage model (Mazars et al., 2015) that has non-differentiable operators in themodel. The first solution is to assume a radial loading path, the second is to analytically regularize the non-differentiable functions and finally is the use of numerical differentiation. The numerical implementation isperformed in Code_Aster finite element software. The performance of the proposed approaches are thencompared in terms of convergence rate, robustness, and computational time at both the elementary level andstructural level. In elementary level, all methods present quadratic convergence for radial loading, but onlythe numerical estimation keep this convergence rate for any loading trajectories. In structural level of areinforced concrete beam in bending with local model, the occurrence of micro snap-backs due to crackingof an element introduces convergence difficulties and hence further strategies are needed.