Influence of Gradual Damage on the Structural Dynamic Behaviour of Composite Rotors: Simulation Assessment

Fibre-reinforced composite structures under complex loads exhibit gradual damage behaviour with degradation of effective mechanical properties and change of their structural dynamic behaviour. In case of composite rotors, this can lead to catastrophic failure if an eigenfrequency is met by the rotational speed. The description and simulation analysis of the gradual damage behaviour of composite rotors therefore provides the fundamentals for a first understanding of complex and partially-unpredicted structural phenomena. Therefore, a simulation tool is developed using a finite element model, which calculates the damage-dependent structural dynamic behaviour of selected composite rotors considering both damage initiation and in-plane damage evolution due to a combination of out-of-plane and in-plane loads. Damage initiation is determined using failure criteria, whereas the gradual damage evolution using a validated continuum damage mechanics model. Numerical results are compared with experimental results for rotor-typical stress states to assess the model quality, which could be later used for damage identification approaches.


Introduction
The development of innovative and efficient damage identification methods is a complex engineering challenge. It requires an interdisciplinary approach that uses experimental, simulation and information-based methods, which encompass the material and structural behaviour. In this context, the understanding of the gradual damage behaviour and of the corresponding structural vibration response of composite structures plays a central role.
A general definition of damage in composite structures as it relates to damage identification methods is given by SOHN [1] as changes in the geometry or in the material properties of a composite structure that affect its current or future performance. In this context, in order to accurately define damage, a comparison between two different structural states has to be performed, one of which is assumed to represent the initial and undamaged state. Therefore, a physical and reliable simulation of the composite damage and the occurring structural phenomena is an important prerequisite for a successful damage identification.

Motivation
An important issue within damage identification is the complexity of composite structures combined with a vast amount of possible combinations of failure modes, damage extents and their positions. Therefore, a common approach for the development of damage identification methods is the application of machine learning algorithms, in which these algorithms are trained to learn datasets containing representative damage configurations [2]. However, in most instances the generation of such learning sets is not realistic in an experimental way due to prohibitive time and cost factors [3,4]. Therefore, a solution to this problem is the development of reliable simulation models for the generation of learning sets considering only a minimum amount of experimental data for fitting purposes [5]. In this context, phenomenological-based damage models capable of predicting different diffuse and discrete damage phenomena have been further developed to describe the progressive damage of composite materials until final failure [6].

State-of-the-Art
The stress-strain behaviour of fibre-reinforced composites can be modelled by constitutive relations based on the generalized Hooke's law in the area of small displacements and under the limits of damage initiation, which is normally the first-ply inter-fibre failure [7]. For the initial failure of fibre-reinforced composites under static, cyclic and highly-dynamic stress, descriptive damage models exist based on the failure criteria of HASHIN, PUCK and CUNTZE [8][9][10]. Furthermore, for the theoretical description of the gradual damage behaviour of composite materials, these advanced models have been further developed and confirmed. The description of damage as a continuous process forms the basis for the further analysis of the gradual damage behaviour of composites and the corresponding damage-dependent dynamic behaviour.

Gradual Damage Behaviour of Composite Materials
For the description of the gradual damage behaviour of composites, caused by operational loads or unpredicted loads such as impacts, verified damage mechanics models for composite materials are already available [11]. Typical failure modes for composite rotors are mainly inter-fibre failure from in-plane loads as has been shown in [12] and a mixture of inter-fibre failure, delaminations and fibre failure from unexpected impact loads [13]. However, the application of these damage mechanics models for rotor-typical loading conditions under consideration of the gradual damage behaviour and the resulting structural dynamic behaviour has not been thoroughly investigated.
An example of such an approach can be a phenomenological damage mechanics model [6] combined with a simulation model [14]. Based on such a model, the gradual damage analysis and the resulting dynamic behaviour of composite rotors can be described. In particular, some first qualitative results have been presented [12,15] for typical fibre architectures which have previously been extensively investigated with their material properties characterised in [16,17].

Damage Initiation Modes
For the calculation of the damage initiation under complex loading conditions, a reliable damage mechanics model was applied that considers the material symmetries by the application of invariants [10,18]. A unidirectional lamina was homogenised to a "material" and separated strength criteria are allocated to each failure mode and to each associated basic strength. The Failure Mode Concept (FMC) focuses on two expressions of the theoretical prediction of failure in composites. The first expression is the derivation of failure conditions for a unidirectional lamina with the prediction of initial failure of the embedded lamina. The second expression is the treatment of a non-linear, progressive failure of three-dimensionally stressed laminates until final failure [10].
According to the FMC, five different failure modes can be identified, as shown in Figure 1. At low load levels, there is formation of three different matrix failure types, IFF1, IFF2 and IFF3, under tensile stress σ 2 > 0, shear stress τ 21 and compressive stress σ 2 < 0, respectively. Increasing stresses result in extended matrix failure and in the formation of delaminations. With further increase of tensile stress σ 1 > 0 and compressive stress σ 1 < 0, resulting in extensive distributed damage, there is the formation of tensile FF1 or compressive FF2 fibre failures, respectively [18].   For every single failure mode, an equivalent stress applying the FMC exists with σ * equ with * = σ, τ, ⊥σ, ⊥τ, ⊥ based on formulated invariants and the strengths from the defined material efforts E f f * . The used formulations for the initiation of damage as well as the final failure are based on the following equations of single efforts [10]: The final effort E f f tot results from a serial connection of the single efforts using the rounding coefficient m [18].

Damage Evolution Approach
For the description of the non-linear gradual damage behaviour, damage evolution laws are then applied, which are based on the constitutive equations for a damaged layer using the applied continuum damage mechanics model. Damage evolution laws have been developed based on FMC using a damage tensor D, which is chosen according to identified degradation mechanisms [6].
A typical example of the gradual damage behaviour of a [0 • /90 • //90 • /0 • ] glass fibre-reinforced epoxy is shown in Figure 2. The specimen undergoes a homogeneous tensile loading, which results in a gradual material degradation. The stress-strain behaviour of the specimen can be viewed with regard to the D 11 damage parameter affecting the E 1 stiffness degradation in the longitudinal fibre direction. Five damage states are also identified for the specific composite configuration, illustrating the gradual damage behaviour of composites [19].
For the description of the gradual damage behaviour, fracture mode-specific evolution equations are predestined, which capture both the stiffness degradation and coupled degradations from the direction-dependent stiffness. These relationships were investigated and described, in particular, for textile-reinforced glass epoxy resins [20]. Based on these findings, further approaches were developed with regard to their implementation into the finite element method [21].

Delamination Mechanisms
In addition to the failure and damage processes mentioned above, which are typically described by continuum mechanics models, delamination is the main out-of-plane failure mode. A variety of methods have been developed for the analysis of the damage behaviour of composites under impact loads causing delamination until final failure [22]. Such methods cover a detailed representation of delamination using cohesive models or the virtual crack closure-integral method, in particular for Mode I and Mode II delaminations [23,24]. Mode III is often not modelled in composites, as it rarely occurs and it is experimentally difficult to characterise. Furthermore, simulation techniques for multiple failure modes, including delamination, have been successfully applied in a simpler, yet reliable method for multiple composite structures [25]. These methods follow material degradation rules by the simulation of the progressive damage propagation at composite structures.
Commonly, the delamination process in composites is described in terms of delamination initiation and delamination propagation, caused by the out-of-plane tensile stress σ 3 and shear stresses τ 13 , τ 23 [26]. Specifically, two different mechanisms of delamination initiation exist, namely, tensile delamination induced by out-of-plane tensile stress σ 3 , and shear delamination under τ 13 , τ 23 , that can initiate the delaminations and their propagation, as shown in Figure 3 for three delamination modes.

Mode I (opening)
Mode II (sliding shear) Mode III (tearing shear) Delamination initiation is mainly treated within the framework of continuum damage mechanics approaches, using quadratic delamination criteria, as shown in Equation (3), where R t 3 is the tensile interlaminar normal strength, R c 3 the compressive interlaminar normal strength, R s 13 the interlaminar shear strength for τ 13 shear stresses and R s 23 the interlaminar shear strength for τ 23 shear stresses [27][28][29].
Different methods have been developed to model delamination within a finite element analysis [22,30]. Firstly, for Mode I and Mode II delaminations, a detailed representation of delamination can be achieved using cohesive models or the virtual crack closure-integral method [23,24]. Mode III is often not modelled in composites, as it rarely occurs and it is experimentally difficult to characterise. Secondly, to model the delamination in a simple, yet reliable method, a different approach has been adopted for several composite structures [25,31]. This approach follows specific structural degradation rules resulting in the degradation of the effective mechanical properties in the delaminated area, and accordingly, the values of Q ij [30,32].

Aim and Outline of the Paper
The aim of the paper is to develop a parametric simulation model capable of depicting the structural dynamic behaviour of composite rotors under propagating damage due to a combination of out-of-plane and in-plane loads.
An elementary rotor geometry is selected with a representative fibre architecture for a CARTESIAN-orthotropic material behaviour. Based on experimental results and their physically-based interpretations [33], a finite element model is developed, which is used for a parametric simulationbased study of the damage-dependent structural dynamic behaviour. The numerically obtained vibration responses are compared with experimental results for rotor-typical stress states.

Selection of Representative Fibre Architecture for Composite Rotors and Manufacturing
For the experimental investigation of the damage behaviour of composite rotors and their resultant dynamic behaviour, a characterisation of failure-tolerant fibre-reinforced composites is performed. In particular, a multi-ply multi-axial fabric is selected resulting in an in-plane CARTESIANorthotropic behaviour.

Multi-Ply Multi-Axial Fabric
The selected fibre architecture is composed of a glass-fibre, non-crimp, multi-ply and multi-axial fabric (NCF) [34,35]. The fabric reinforcement has an area density of 1.90 kg/m 2 and a ply thickness of 1 mm. The composite lay-up consists of four such fabrics [(0 Figure 4, resulting in a total laminate thickness of 4 mm and in an in-plane orthotropic behaviour. The inner and outer diameter of the rotor are 60 mm and 500 mm, respectively. This material was extensively investigated in previous research projects, where the results indicated a gradual damage behaviour [17]. The manufacture of the composite rotors is carried out using the vacuum-assisted resin transfer moulding (VARTM) process [17].

Quality Assessment of Manufactured Rotors
In order to determine the thickness variation, quality assessment is performed on the manufactured NCF-rotors, using the optical marker recognition system TRITOP and the digitizing system Advanced Topometric Optical Sensor (ATOS) from the company GOM. Multiple measurements are taken for each rotor and combined using the system TRITOP, in order to minimise the error by the data combination. Each single 3D scanning of the entire test object has a spatial point-to-point distance of approximately 0.27 mm by using a camera system with an optical resolution of 1280 × 1024 pixel, enabling the measurement of small, local deformations and manufacturing defects at the rotor surface.
The rotor is placed on a rotary table and it is measured with a high resolution of approximately 30 µm, as shown in Figure 5. For the determination of the rotors thickness variation, the acquired experimental data are further analysed using parametric fitting, in order to determine the shape of the rotors. A quadratic polynomial curve is used as fitting equation of the type Specifically, 36 evenly distributed rotor sections are selected with approximately 650 data points each. Then, the mean thickness value and the standard deviation is calculated along the radius of the rotor for 20 evaluation points, and a quadratic polynomial curve is fitted for each manufactured rotor. The resulting manufacturing variations of the rotors caused by the internal pressure are presented in Figure 6.

Development of a Parametric Simulation Model
A simulation model is developed, which predicts the damage-dependent structural dynamic behaviour of the investigated composite rotors under consideration of both the damage initiation and the in-plane damage evolution from a combination of out-of-plane and in-plane loads. Within this model, the damage initiation is determined by means of the failure criterion. The gradual damage evolution can be described with a validated continuum damage mechanics model, which captures the degradation processes after considering the failure modes and respective damage coefficients.
In addition, the simulation results are validated with experimental results from intact and damaged rotors regarding the gradual damage behaviour and the corresponding damage-dependent dynamic response [33].
For the development of a parametric simulation model, the finite element (FE) software ABAQUS is used, and a parametric PYTHON script is implemented for the investigation of the damage-dependent dynamic behaviour. This enables the generation of a FE model with parameters that define the position and magnitude of the out-of-plane load. Subsequently, a finite element analysis (FEA) is performed, and a detailed investigation of the gradual damage evolution as well as the resulting dynamic behaviour of the rotors is conducted. The flow chart of the developed simulation approach is shown in Figure 7.
For each type of rotor, the geometry is defined according to the measured averaged thickness of the manufactured rotors, which was previously experimentally determined. Specifically, the preform thickness and the internal pressure-induced deflection of the mould are taken into account. Manufacturing imperfections such as the deviation of the fibre orientation or local resin pockets are not modelled. Material parameters determined in previous investigations are also retrieved [17]. The model is meshed using continuum three-dimensional, 8-node linear brick solid elements with incompatible modes of type C3D8I with approximately a total number of 15,500 nodes and 7500 elements as shown in Figure 8. Based on the selected type of element, results from previous investigations as well as the comparison of the numerically estimated to the experimentally measured modal properties, one element in the thickness direction is sufficient for this kind of analysis and the desired resolution. The quality of the mesh plays a very important role both at the damage prediction as well as to the eigenfrequency calculations. In order to identify an optimal size of the elements that result in sufficient good results, a convergence investigation for the first eigenfrequencies has been performed at preliminary investigations. This investigation showed that with this element size we achieve a convergence. A qualitative comparison between experimentally and numerically determined damage states took place. A further more detailed investigation for the mesh quality will be taken into consideration in future investigations for similar rotors with a variable thickness along the radius. For the boundary conditions of the model, the degrees of freedom located at the inner area of the rotor are constrained in order to simulate the clamping area. In the boundary nodes all the degrees of freedom are constrained. The clamping is modelled by selecting the nodes located at the inner area of the rotor, and then by setting the translational degrees of freedom of these nodes in the out-of-plane direction of the rotor as well as the respective rotational degrees of freedom to zero, as shown in Figure 8.  A user-defined field subroutine (USDFLD) [21] is included to model the damage initiation and the gradual damage behaviour of NCF composite rotors. This subroutine enables the description of the damage initiation and propagation, with respect to previously developed phenomenological damage models [20]. In this work, the modelling of delamination using a reduction of effective mechanical properties is followed with regard to Equation (3). This technique, although simple, is considered as sufficient, according to experimentally gathered information from ultrasonic testing and its comparison to the simulation results.
The considered loads are derived from previous experimental investigations [33]. The out-of-plane load, which creates the initial damage, is modelled as a distributed load at the nodes included in the contact area between impactor and rotor [33]. The rotor run-up is modelled as rotational velocity steps of 1000 rpm in a range from 8000 rpm through 14,500 rpm for a NCF rotor based on previous experiments [33]. The applied increasing body forces cause the increase of complex states of stress as well as the initiation and the resulting damage evolution.
In order to consider the manufacturing-relevant residual stresses of the rotor, the curing process is simulated as a pre-stress effect, where a thermal static analysis is included with a temperature change of ∆T = 60 K. Subsequently, the stress fields are calculated and used for the failure analysis as well as for the calculation of the damage initiation. Based on the calculated data, the damage evolution at each layer is estimated and the corresponding effective material properties are included in the model. The estimated developed damage is subsequently transferred as starting state for the next loading step.
For each estimated damage state, the damage-dependent dynamic response of the rotor is calculated using a numerical modal analysis and, approximately, the first 20 eigenmodes are estimated for a frequency range from 20 Hz through 1000 Hz. Subsequently, a random response analysis is performed for the calculation of the power spectral densities at selected degrees of freedom [33]. The excitation is simulated as an uncorrelated band-limited white noise in the frequency range of interest.
The consideration of damping in commercially available software programs is state of the art for classical isotropic monolithic materials. However, the direction-dependent damping coefficients of anisotropic materials are not supported. To overcome this deficiency, a post-processing feature has been developed in MATLAB, able to calculate the modal loss factor of anisotropic fibre-reinforced composite structures and is here also implemented [36,37]. The model includes the anisotropic material damping by an extension of the FEM using the strain energy concept [38,39]. The dynamic response is extracted at 128 degrees of freedom corresponding to the measurement points of previous experimental modal analysis [33].
The respective damage-dependent change of the material damping is not considered due to the missing material damping values. The estimation of the damage-and direction-dependent damping coefficients of anisotropic materials is still a challenge, which is already in the scientific focus [40], but still has to be addressed in future investigations.
The described simulation analysis is performed sequentially for each defined rotational velocity. The simulation delivers numerical results for the local stresses of the rotor, for the resulting efforts and the respective failure modes as well as for the dynamic response of the rotor.

Assessment of the Parametric Simulation Model
A comparison between the experimental and the numerical results is performed for the validation of the simulation approach [33]. First, the simulation model is assessed with regard to the damage state and the damage evolution. Then, the dynamic vibration response is compared and the modal dynamic properties, i.e. the mode shapes, the eigenfrequencies and the modal loss factor are compared for multiple structural states. Finally, the dynamic behaviour is considered under the effect of the evolving damage and a comparison between experimental and numerical change of eigenfrequencies is conducted.

Damage State and Damage Evolution
A comparison of the damage effect of an out-of-plane compression load of 16 kN an area between experimental and numerical results reveals a similar damage extent as well as the same type of damage, as shown in Figure 9. The contact area was simulated with a diameter of 10 mm, based on experimental results [33]. The same main failure modes appear both in the model and in the experiments, which are inter-fibre failure, fibre-matrix shear failure and the formation of delamination. Overall, the investigation of the damage state at different applied compression loads confirms a good qualitative agreement of the damage extent according to experiments and simulation.

Delamination
Inter-fibre failure The damage initiation and propagation due to an in-plane load under a rotational velocity is shown in Figure 10. The damage evolution is not-symmetric due to the selected fibre architecture and the resulting in-plane CARTESIAN orthotropic material behaviour. A qualitative accordance is observed with regard to the damage extension as well as the main failure modes that are apparent, which are inter-fibre failure due to tensile and shear stresses.

Modal Properties
With regard to previous experiments [33], the developed model provides all the necessary information for the dynamic behaviour of the investigated composite rotors, mainly the mode shapes, the eigenfrequencies, the consequent change of eigenfrequencies due to damage increase, the modal loss factor and the single as well as the mean power spectral density. For the model validation, the experimental and numerical properties are compared for a frequency range from 20 Hz through 2000 Hz, including 16 experimentally measured eigenmodes.

Mode Shapes of Intact Rotors
A qualitative comparison between a calculated and an experimentally determined mode shape (3,1) for an intact NCF-rotor is shown in Figure 11. The calculated mode shape corresponds with the experimentally measured one, providing a first qualitative verification of the simulation model. For the quantitative evaluation of the developed model with respect to the mode shapes, the Modal Assurance Criterion (MAC) is applied, which yields a good statistic indicator and a degree of consistency between mode shapes. The MAC is calculated as which forms the linear dependency of two modal vectors x exp , x num . Each modal vector x exp , x num is the vector of a mode shape with a length of n = 128 measured and calculated points, respectively. A MAC value of unity represents a collinearity of the two vectors, which indicates fully consistent mode shapes. In contrast, a MAC value of zero indicates the orthogonality of the two eigenvectors, which means that the modes are not consistent. A tolerance limit of 70% is usually acceptable for a consistent correlation between the investigated mode shapes [41]. A high agreement of the mode shapes is observed for the 11 identified EFs for the NCF I.3 rotor, as it is shown in Figure 12, as well as for the other investigated rotors.

Modal Loss Factor
The comparison of the numerical and experimental results of the modal loss factor is shown in Figure 13. The quantitative comparison showed a good accordance with deviations up to 4%. The differences between experimental and calculated modal loss factors are due to the contribution of friction-caused damping at the clamping area.  A mean modal loss factor of 1.61% and 1.12% is experimentally and numerically determined. Due to the small changes, a constant modal loss factor of 1.12% will be used for the simulation. For more complicated structures where the mode shapes are substantially different, a variance of the modal loss factor is expected and in this case a vector array should be used.
For the damaged states, the modal loss factor remains almost constant, with a mean value of 1.74% with a standard deviation of 0.33%, based on experimental data from all investigated rotors. These values are similar to the modal loss factor determined previously for the undamaged states. Therefore, the modal loss factor from the undamaged state can be used as constant values for both damaged and undamaged states, although it is expected that the damaged states should have an increased modal loss factor due to friction in the interface of the inter-fibre cracks. The possible friction in the interface of the inter-fibre cracks can influence the damping behaviour of the investigated composite rotors. The selected method of extracting the modal damping ratio provides a first assessment of the global frequency-dependent damping [42]. When a further more detailed damping estimation is required, then an estimation of the spatial distribution of the damping properties using a dynamic mechanical analysis could be applied [40]. In this way also damage-caused increase of the anisotropic material damping can be investigated, which originates from dry friction in fibre-reinforced polymer-based composites [43]. Figure 14 shows the experimentally and numerically determined eigenfrequencies of a NCF rotor and two eigenfrequencies, the (2,0) and (2,1) are exemplary presented. Specifically, one with a non-monotonic change of the eigenfrequency (left) and one with a monotonic decrease of the eigenfrequency (right) is observed.

Change of Eigenfrequencies Due to Damage Increase
The relation between applied in-plane loads and the corresponding shift of eigenfrequencies due to damage increase is investigated for five different damage sequences. For each damage sequence, a total of 16 eigenfrequencies are investigated. For each type of rotor, a normalised in-plane load F c is introduced, where each rotational velocity n rpm is divided from the maximum achieved rotational velocity n max rpm for each sequence.   Figure 14. Experimentally and numerically determined change of the eigenfrequency for two eigenmodes of a rotor, with and without initial damage.
A quantitative assessment of the model's performance is conducted with regard to eigenfrequency change. The correlation is selected as a statistical criterion for this assessment. The correlation is calculated between the mean fitted curves of the experimental (EF exp ) and numerical (EF num ) eigenfrequencies, and the results are shown in Figure 15. There is approximately more than 80% correlation between experiment and simulation for all investigated eigenfrequencies.

Conclusions
Based on experimental results and their physically-based interpretation [33], a parametric simulation model is developed and verified, which numerically determines the damage-dependent structural dynamic behaviour of the investigated composite rotors. This model can be later used, e.g., for the generation of input data required for different damage identification methods. A very good agreement between measurements and simulation results is identified and validated for rotor-typical stress states. Overall, the good compliance of experimental measurements and numerical simulations shows that the proposed numerical approach can provide a reliable way for simulating the dynamic behaviour of composite rotors with a complex fibre architecture under different structural conditions.
The deviation between experimentally and numerically determined results is reasoned by production-related features, such as the thickness variation of the manufactured rotors, preform thickness variations, small changes to fibre orientation and local resin pockets. On the one hand, a better manufacturing process would result in a higher reliability and therefore in a better accuracy of the models. On the other hand, an exemplar-specific modelling and simulation approach of the rotors could also improve the accuracy of the models and could be applied in the field of digital twins of composite rotors. Funding: The authors would like to express their gratitude for the financial support from the Deutsche Forschungsgemeinschaft (funding code GU 614/14-1).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: