In a process of permanent innovation and increase of turbojet efficiency, the low pressure turbine is designed to run at higher rotational speed in future engines. Such an evolution requires a precise evaluation of their effects on dynamic behavior of the turbine. Low-pressure turbine blades in aeronautical industry are generally composed of a shroud, whose main purpose is to keep the tightness of the primal air stream, but it produces in consequence a mechanical damping into the dynamic response of the full blade turbine assembly through frictional contacts. Depending on the industrial design choices, the geometry of these blades can induce a pre-torsion effect which will straighten the shrouds during the assembly process, in this way they are brought into contact and the blades are pre-loaded. Predicting the dynamic behavior of such structure and calculating forced responses is a major challenge considering the complexity involved. The numerical model used in this article is based on the geometry of a test bench which was designed to reproduce the pre-torsion assembly loading on the blades through a geometrical interpenetration, and to evaluate the dynamics of the contact between the shrouds. It includes only three blades with a simplified geometry but represents the main dynamic characteristics of a real turbine. A static analysis on a commercial FE software of the test bench blades assembly has been carried out to evaluate the contact area on the shrouds, and determine the deformation related to the pre-torsion loading. Modal analyses were also realized to determine the natural frequencies of the structure and the modal shapes according to the imposed contact status. Then a reduced order model of the structure is generated, and is used into a specific numerical code. A procedure is implemented to handle the interpenetration generated by the pre-torsion of the blades, which is the key part of the static loading of the structure. In order to solve this non-linear dynamic problem in the frequency domain, based on the static loading results, the Harmonic Balance Method (HBM) is coupled to the Dynamic Lagrangian Frequency Time (DLFT) method or a classic penalty method to evaluate the non-linear contact forces. Numerical results around the first and third bending modes of the structure are computed and a parametric study is presented, pointing out the impact of the static normal load distribution depending of the method used, the friction coefficient between, the excitation force level, etc. The precision of the nodal coupling scheme used to manage the contact between the shrouds meshes for the numerical simulation is also evaluated by changing the meshing, in order to quantify such hypothesis.