The mechanisms for nonlinear saturation of a bluff-body stabilised turbulent premixed flame are investigated using LES with the transported flame surface density (TFSD) approach to combustion modelling. The numerical simulation is based on a previous detailed experimental investigation. Results for both the unforced non-reacting and reacting flows are validated against experiment, demonstrating that the fundamental flow features and predicted flame structure are well captured. Key terms in the FSD transport equation are then used to describe the generation and destruction of flame surface area for the unforced reacting flow. In order to investigate the non-linear response of the unsteady heat release rate to acoustic forcing, four harmonically forced flames are considered having the same forcing frequency (160 Hz) but different amplitudes of 10 %, 25 %, 50 % and 64 % of the mean inlet velocity. The flame response is characterised using the Flame Describing Function (FDF). Accurate prediction of the FDF is obtained using the current approach. The computed forced flame structure matches well with the experiment, where effects of shear layer rollup and growth of the vortices on the flame can be clearly observed. Transition to nonlinearity is also observed in the computed FDF. The mechanisms leading to the saturation of the flame response in the higher amplitude case are characterised by inspecting the terms in the FSD transport equation at conditions when the integrated heat release is at its maximum and minimum, respectively. Pinch-off and flame rollup can be seen in snapshots taken at the phase angle of maximum integrated heat release. Conversely, intense vortex shedding and flame-sheet collapse around the shear-layer, as well as small-scale destruction of flame elements in the wake, can be seen in snapshots taken at the phase angle of minimum integrated heat release. The pivotal role of FSD destruction on nonlinear saturation of the flame response is confirmed through the analysis of phase-averaged terms in the FSD transport equation taken at different locations. The phase-averaged subgrid curvature term is found to concentrate in the cusps and downstream regions where flame annihilation is dominant.