The aim of the present study is to analyze the saturation regime of the Toroidal Alfven Eigenmodes (TAE) in the LHD plasma, particularly the MHD burst. The linear and non linear evolution of the TAEs are simulated by the FAR3d code that uses a reduced MHD model for the thermal plasma coupled with a gyrofluid model for the energetic particles (EP) species. The linear simulations indicate the overlapping of $1/2-1/1$, $2/3-2/4$ and $3/5-3/6$ TAEs in the inner-middle plasma region and frequency range of $45-75$ kHz, triggered by EPs with an energy of $T_{f} = 45$ keV and EP $\beta = 0.022 $. The non linear simulations show that $2/3-2/4$ and $3/4-3/5$ TAEs are further destabilized due to the energy transfer from $1/1-1/2$ TAE, leading to a broad TAEs radial overlapping and the MHD burst triggering. The energy of $1/1-1/2$ TAE is also non linearly transferred to the thermal plasma destabilizing the $0/0$ and $0/1$ modes, inducing the generation of shear flows and zonal currents as well as large deformations in the thermal pressure and EP density radial profiles. The non linear simulation reproduces the same succession of instabilities and the same frequency range with respect to the experiment. The instability propagates outward during the bursting phase, showing a large decrease of the EP density profile between the middle-outer plasma, pointing out the loss of part of the EP population that explains the decrease of the plasma heating efficiency observed during the MHD burst.