Computational modelling of dynamic delamination in morphing composite blades and wings

Morphing blades have been promising in lifting restrictions on rated capacity of wind turbines and improving lift-to-drag ratio for aircraft wings at higher operational angles of attack. The present study focuses on one aspect of the response of morphing blades viz. dynamic delamination. A numerical study of delamination in morphing composite blades is conducted. Both components i.e. the composite part and the stiffener are studied. The eXtended Finite Element Method (XFEM) and nonlocal continuum mechanics (peridynamics) have both been used to study fracture in the isotropic stiffener used in conjunction with the blade. As for the composite morphing blade, cohesive elements are used to represent the interlaminar weak zone and delamination has been studied under dynamic pulse loads. Intraply damage is studied using the nonlocal model as the peridynamic model is capable of addressing the problem adequately for the necessary level of sophistication. The differences and similarities between delamination patterns for impulsive, dynamic, and quasi-static loadings are appreciated and in each case detailed analyses of delamination patterns are presented. The dependence of delamination pattern on loading regime is established, however; further parametric studies are not included as they lie beyond the scope of the study. Through the use of fracture energy alone the nonlocal model is capable of capturing intraand interlaminar fractures. The proposed modelling scheme can thus have a major impact in design applications where dynamic pulse and impact loads of all natures (accidental, extreme, service, etc.) are to be considered and may therefore be utilised in design of lightweight morphing blades and wings where delamination failure mode is an issue.


INTRODUCTION
As the level of functional requirements for an aircraft or turbine blade elevates so does the level of sophistication in design. One such scenario emerges when higher efficiency in energy production or higher pitch controllability to increase lift, minimise drag, and delay stall are desirable. To devise a conceptual design for a realistically functional blade one must consider several requirements. Not only must the blade be aerodynamically efficient, among other things, but also structurally safe and operative. As such, the material for the blade must possess high strength-to-weight and stiffness-to-weight ratios but the design must also be flexible enough to render possible controllability thus a balance between competing goals must be reached. On the one hand, the structure must be robust and reliable under the application of possible loads, on the other hand, the shape of the blade must be aerodynamically optimised to maximise lift-to-drag ratio and avoid stall (thus loss of lift) at operative angles.
The primary structural requirements for the blade are its safety when subjected to service and ultimate loads during its lifetime as well as its fatigue durability [1][2][3]. To ensure structural robustness, optimal material and thickness distributions along the axis of the blade must be considered. Given the state of knowledge has grown in parallel with computational power, the designs of blades have changed throughout their development history [3,4].
Advances in structural aerodynamics have given rise to constant evolution in novel blade designs of improved airfoil cross sections and wings with continuous shape morphing. This allows these flexible structures to adapt continuously their shapes with respect to aerodynamic conditions and applied loads (see Fig.1). The smooth shape changes of the airfoil feature a strong potential to increase the lift capability and the aerodynamic stability (see Fig.'s 2 and  3). This morphing concept mimics the biologic systems in nature and currently attracts the attentions of researchers from both academia and industry. Hence, morphing airfoils have become a valuable candidate for innovative development in both aircraft and wind turbine technologies. Blades are nowadays designed employing either the active control concept (using solid state transducers, electro-active polymer, piezoelectric actuators) or passive control concepts. More specifically, some new wind turbines are considered with passive blades shape adaptation (without actuators, sensors, or any control devices). The simplicity of this passive pitching control system makes it a low-cost addition to the device and could dramatically increase the energy conversion of the rotor. As such, several recent works have been conducted on the aerodynamic characterization of such blades and experimental and computational studies galore have contributed to the field [3,[5][6][7]. However, the mechanical degradations associated with operations in designing the shape of morphing for the blade and structural limitations which arise have not been completely addressed.    [8] Morphing wings and blades have been proposed by researchers ever since the dawn of aviation when the idea of morphing an aircraft's wing characteristics through continuous, rather than discrete, movable aerodynamic surfaces was believed to hold the promise of more efficient flight control [9]. Morphing technology has been proposed not only for large-scaled manned aircraft, but also for Unmanned Aerial Vehicles (UAV's) where morphing combined with advances in actuator/sensor and material technology has given rise to a surge in interest in alternative morphing configurations capable of matching multiple mission profiles encountered sequentially in a single flight mission [10,11].
Composites have been proposed and used for the sake of aeroelastic tailoring [12]. As composites possess the potential to be stiffness-tailored several considerations must be made in their use. A multilevel approach accounting for aero-structure interaction might be needed to acquire an understanding of the actual spatial and temporal distribution of loading on the wing/blade. As such fluid-structure interaction (FSI) might be needed to be taken into account for accurate derivation of loading. This has been done by researchers [13]. These analyses, usually conducted using an arbitrary Eulerian-Lagrangian (AEL) finite element method, can quantify service and accidental extreme loads if enough information is provided. Such analyses have been done for metal and composite structures successfully [14,15].
While, more often than not, these blades and wings are subjected to service loads it may happen that ultimate loads of static or dynamic nature may be applied to them. As such, some or all modes of failure may be triggered. These include fibre failure, matrix cracking, in-plane shear failure, fibre micro-buckling, fibre-matrix debonding, and delamination. While most modes of failure occur within a ply of composite (intralaminar modes) delamination transpires in between plies and is thus referred to as interlaminar mode of failure. The resin rich layer (RRL) between adjacent plies where no fibres are present provides a weak link for delamination to initiate and propagate.
Delamination in composites due to dynamic loads has been widely studied [16][17][18]. One aspect of delamination in pulse-loaded laminated plates is attributed to 'early time' response i.e. spalling or scabbing failure modes, which transpire within micro-seconds, where the through-thickness compressive waves give rise to tensile stress waves on reflection from the free back face and result in mode I interlaminar failure [19][20][21][22]. Plate action which follows could lead to mode II and mixed-mode interlaminar fracture [23,24] similar to dynamic debonding of soffit plates in a strengthened beam [20].
Detailed analyses of such cases using numerical methods, such as the finite element method, and analytical models, as first-and higher-order shear deformation laminated plate theories [25,26], are possible and the theories can be used to study such cases. Analytical models may lead to closed-form solutions only for simple geometries, idealised material and loading conditions, and specific boundary conditions, and are thus limited. Furthermore, while initiation of cracking in the interlaminar zone is accurately predicted in these cases, crack propagation is a nonlinear problem and poses numerous new challenges which may only be adequately addressed within a numerical framework. Numerical modelling schemes as finite element method lead to a comprehensive framework for simulation of such cases as they allow for incremental solution of the nonlinear problem and lead to accurate results. Alternative powerful methods as nonlocal continuum mechanics-based model of peridynamics have also proven capable to capturing crack propagation, branching and coalescing accurately. The aim of the current study is thus to ensure a framework exists based on these methods which can address adequately all the major meso-scale and macro-scale phenomena as far as different modes of fracture are concerned. Modes I and II, and mixed mode fracture are included and the primary hypotheses of each model explained. Interlaminar and intralaminar fracture are considered and limitations discussed so the framework is capable of addressing the problem universally. The outline of this work is as follows.
In section 2 several analyses are conducted on both the stiffener and the morphing blade. Static analyses and frequency linear perturbation method are used to study the response for each case as well as to extract eigenvalues. Extended finite element method (XFEM) and nonlocal continuum mechanics model (peridynamics) are used to look into the maximum load carrying capacity of the stiffener and a linear perturbation eigenvalue solver is utilised to extract natural frequencies and modes. In section 3 dynamic delamination using cohesive traction-separation law and nonlocal continuum mechanics (peridynamics) have been used to study delamination in the composite morphing blade under pulse loading. Several regimes of loading are considered, however, intra-ply damage has not been included as (1) the focus has been on delamination and it is believed that this mode of damage precedes fibre failure, which is the detrimental mode of damage to a composite lamina and (2) enrichment functions (asymptotic crack tip functions) are not implemented for anisotropic media and the use of continuum damage mechanics requires implementation of a mesh-objective damagemechanics based constitutive model which lies outside the scope of the present study. The intraply damage mode may be included using peridynamics as shown in this work. Section 4 concludes the studies conducted.

ANALYSES
Airfoil cross sections in conjunction with twisted-pitch blades have been traditionally used as the preferred geometry in the design of wind turbine blades and aircraft wing structures. These forms have an aerodynamically advantageous performance as they generate a high ratio of lift to drag. Composite materials along with wood are used in the construction of blades as they possess superior strength-to-weight and stiffness-to-weight ratios. Glass fibre reinforced plastics (GFRP's) also have very good fatigue resistance important to blade root design. As for flexible wings, carbon fibre reinforced plastics (CFRP's) and aluminium are used as materials of preferred characteristics. As for geometries of the cross sections usually aerodynamically optimised sections are adopted from the libraries of NREL for wind turbine blades, and NACA for aircraft wings.
The stiffeners are usually sections without holes (see Fig. 4a) while the blade is made of a hollow section to allow morphing flexibility (see Fig. 4b). Both models, however, contain a hole to allow them be connected to the spar bar of the main wing structure.
As far as the global structural response is concerned both finite element simulations and analytical models may be used to acquire an understanding of static and dynamic behaviours of the blade. The level of complexity for analytical models depend on the desirable level of details to be derived accurately and economically which depends on what the intended use of analytical model is and the features it needs to represent faithfully. Simplified methods of analysis allow the derivation of some characteristic response parameter (e.g. displacement at a particular point) using a generalised parameter model such as an SDoF model [27] and/or more complicated models [28,29]. As an example the entire morphing blade or stiffener may be modelled as a non-prismatic beam semi-rigidly connected to a support as depicted schematically in Fig. 5a, a model which in itself may be viewed as the superposition of two models of Fig.'s 5b and 5c in static analysis or as a simplified generalised parameter SDoF model to be analysed dynamically.
Despite the fact that this model provides information on global response it does not provide the details of response at each section level nor does it allow for nonlinear damage phenomena as fracture, damage, or delamination.
As far as aerodynamics is concerned, the lift increases as a function of the angle of attack up to a limit, where loss of lift is encountered (see Fig. 6, which gives the dimensionless lift coefficient CL for a rigid airfoil as a function of the angle of attack). From the aerodynamics point of view the morphing blade is shaped in such a way that it can increase the angle of stall through flexible morphing thus to delay stall and provide lift for larger angles relative to air velocity, which is beneficial to energy production (for turbines) and flight mechanism (for aircraft).

Flexible beam
Built-in support (c) Fig. 5: (a) Schematic of a non-uniform beam semi-rigidly connected to a pin support, (b) a rigid beam of the same cross section as (a) and semi-rigidly supported, (c) the flexible non-uniform beam rigidly connected to the support. (Note: superposition must be done as if these were springs in series i.e.  [30]. In the sequel, we shall present at first the results of static analyses on the stiffener, followed by analysis of failure using both extended finite element method (XFEM) [31][32][33][34] and peridynamics [35][36][37][38].

Static analysis
As mentioned above, the stiffener is a non-prismatic beam with no holes connected to a spar bar which is sometimes controllable in length [9]. This section deals with static analysis of the stiffener, made of aluminium, and looks into the load vs deflection when subjected to uniform lower surface pressure load, replicating the operational pressure distribution. The effect of geometric nonlinearity is once excluded and once included in the model and the differences are discussed. In the latter case, the equations of equilibrium (and by extension, the equations of motion in the dynamic case) are derived based on the deformed configuration of the body undergoing deformation as a function of time. The type of aluminium selected is AA7050-T7451 with the properties as follows: The finite element model has been set up assuming the geometry of the section is according to NACA 8412 airfoil and the chord length is 100mm. A 3D geometry has been created in ABAQUS 2018, however, out-of-plane degrees of freedom are restrained as a plane strain condition is presumed and the hole is fully clamped as the effect of partial constrain of the spar bar has not been investigated. As such, this component (stiffener) is studied in isolation. As depicted in Fig. 7 a uniform lower surface pressure has been applied and the displacement at the trailing edge has been recorded. Fig. 8 shows the finite element discretisation of the stiffener. The out-of-plane thickness is taken as 4mm and there are 4 elements through the thickness.     It shows at large deflections (around the aspect ratio of U2/L=32/100 ~32%) the deviation of nonlinear response from the linear response starts and grows as an ascending function of displacement. Fig. 11 shows a comparison of horizontal displacement at the trailing edge. It is concluded that substantial contraction subsequent to some elongation is observed due to finite rotations at the vicinity of the support and the airfoil section folding onto itself. The graphs in Fig. 11 are presented as displacement vs step time which is a parametric form of representation. The reason for this is that force vs displacement in this case would lead to a nonphysical representation and will not be useful.

Frequency analysis
Along with the static analyses conducted, the eigenvalue problem has been solved to obtain the natural frequencies and natural modes of vibration. This is useful for two reasons: (1) it helps identify the regime of pulse loading as impulsive, dynamic, or quasi-static, as this depends on the ratio of loading duration to natural period of the first mode and (2) it may be used to solve for the dynamic response of linear systems using the method of modesuperposition which asserts the response can be obtained as a linear combination of natural modes.
The eigenvalues and eigenvectors may be obtained based on Lanczos or subspace algorithms [40], and the latter has been used here. Fig.13 depicts the first four modes of vibration along with their associated frequencies.
The number of nodes and antinodes in each mode corresponds to the number of that mode and the wave number (inverse of wave length) decreases as the number of the mode (and its frequency) increases. It is expected that the dynamic response is a superposition of weighted modes, however, for dynamic pulses of long duration compared to the natural period of the first mode, most of the response occurs within the first few modes and the inclusion of all modes is not necessary.
Thus, more often than not, the regime of loading is determined by the ratio of pulse duration to the largest natural period i.e. the first mode's natural period. Since fracture introduces a discontinuity into the displacement field and stress/strain field singularities at the crack tip, the analysis of fracture requires special techniques such as continuum damage mechanics (CDM), or extended finite element method (XFEM). Within the classical paradigm of fracture mechanics(such as linear elastic fracture mechanics (LEFM)) singularities are circumvented through the use of energy release rate (G-values or Jintegral) or stress intensity factors (K-values) [41][42][43]. While CDM is based on weakening the failed zone based on the instantaneous amount of energy dissipated, XFEM is based on enrichment of the displacement function using asymptotic crack tip functions and the Heaviside step function [44]. XFEM has been used in the past successfully to study fracture initiation and propagation in isotropic and anisotropic media [31,34,45]. Eq. (2) shows the extrinsic enrichment of the displacement field using the Heaviside step function of Eq. (3a) and the asymptotic crack tip functions for an isotropic medium shown by Eq. (3b).
Where ( ) is the Heaviside function defined as follows: Where is the location of a Gauss point and * the point on the crack closest to , and is the unit outward normal to the crack at * and the asymptotic crack tip functions are as follows: As crack initiation is stress-based and propagation energy based, the stress at failure has been used along with the G-value to replicate the cracking process (table 1).
Under uniform pressure, the response entails the initial elastic response discussed in 2.1.1 followed by initiation and propagation of a discrete crack until ultimate failure as follows: The resistance curve (load-deflection) for uniform vertical traction is shown in Fig. 15. The maximum load in this case is obtained as 527N which occurs when the tip displacement is 3.33mm.
On a closer examination, the progression of the cracking process and the associated stress distribution maybe analysed as follows: At first stress concentration occurs around the square hole with maximum stress occurring at sharp corners. This is depicted in Fig. 16a. Then a crack forms and starts to propagate (Fig.  16b) until the two boundaries are hit by crack tips (Fig. 16c). At this point the structural integrity is lost and failure has occurred.
As level sets are used in XFEM, the contours of − and −levels are given in Fig.'s  17a and 17b. −level sets are defined to distinguish between the areas of enrichment above, along, or below the crack path while -level sets are defined to discern between at, behind, or in front of the crack tip thus > 0 ( < 0) means above (below) the crack path and = 0 means along the crack path, while > 0 ( < 0) means in front of (behind) the crack tip and = 0 means at the crack tip. Evolution of cracking and relaxation of stresses near the cracked zone (shift in the position of the maximum stress due to energy release and drop in stress levels in the vicinity of the crack tip and faces) is governed by energy release rate and is discussed in detail by researchers [46,47]. As reproducing elements contain a crack their status shifts from pristine to failed.  An alternative method which allows for simulation of crack propagation and branching is peridynamics [35]. It is based on a reformulation of nonlocal continuum mechanics and allows for long-range interaction between continuum points. This method is used here to predict fracture patterns in aluminium and in a later section to look into fracture in composite blades and has also been extended to composite morphing blades. At the heart of the theory is inclusion of integrals in the equation of motion of elements, called particles here, which eliminates the need for a continuous displacement field as in the classical continuum mechanics theory [35,48]. In a previous work, the bond-based peridynamic theory was extended to modelling fracture in anisotropic media by using spherical harmonics [49]. The details of the method may be found in the literature [35,50], and to avoid longevity and repetition here we will only state the fundamental equations.
The equation of motion in peridynamics reduces to an integral equation as follows: Where ρ is the mass density in the reference configuration, ( , ) is the displacement vector field, is a pairwise force function, which is of the dimension of the force per volume squared that the particle at position exerts on the particle at position ′ , and is the body force vector field. While it is expected that the integral in Eq. (4) must be evaluated over the entire volume of the medium ( ) it is neither economic nor essential for this to be the case. As such, the concept of a horizon i.e.
= { ′ ∈ , | ′ − | < } is introduced which signifies, as its radius ( ), the finite distance over which two particles can interact and the long-range forces are non-negligible. Defining as = ′ − the relative position of two points in the undeformed configuration and by = ( ′ , ) − ( , ) = ′ − the relative displacement of the two points at time ≥ 0 , the fundamental properties of the pairwise force function are as follows: ( + ) × ( , ) = The deformation in a single bond is defined using the stretch parameter as follows: When subjected to uniform lower surface vertical traction the pattern of fracture will be as in Fig. 19. These simulations are conducted using the in-house code developed at Imperial College London. The maximum load bearing capacity for the model is 516N which is within 2% of the result obtained by XFEM.
It is possible to define a damage index (parameter) as in Eq. (11) similar to that obtained in CDM. However, a detailed study of macro-damage as a result of micro-failures of the bond lies beyond the scope of this work.
The same method will be used in the sequel to study fracture in a composite morphing blade.

Analysis of the morphing blade
The composite morphing blade forms the primary subject of this investigation since the idea of morphing is the core of sustainable flight mechanism or energy production for larger angles of attack. A single ply in a composite laminate comprises fibres and a matrix which may be arranged in several architectures e.g. unidirectional or woven. As such, the ply may be orthotropic or transversally isotropic and the plane of isotropy, in the latter case, depends on fibre architecture. A laminated plate's constitutive matrix is therefore as in Eq. (12): It must be noted that at the material level the general form of a constitutive matrix for an anisotropic material is = . When this relation is integrated through the thickness of the laminated plate once for forces i.e. directly, and then when the moment equation i.e. stress times lever arm is integrated Eq.(12) is obtained.
As the geometry of the blade is complicated, numerical simulations are required to solve for the response under static and dynamic loads and we shall focus on the properties at ply level and will not directly use Eq.(12) to obtain the response. Furthermore, we will only consider unidirectional (UD) composites. Fig. 20 shows the schematic of possible direction of fibres in each ply (direction 1) and transverse direction (2) as well as through-thickness direction (3) for two consecutive plies of the composite laminate. For the sake of this study only carbon fibre reinforced plastic (CFRP=graphite-epoxy) has been considered, however, the method is general and can accommodate other types of FRP's. Table 2 shows material properties used in the analyses. For the sake of comparison, some other composites are also included.
Normally fibre orientations vary through the thickness of a laminate to provide reasonable directional stiffness and strength. For instance, for aeronautical composites, material orientation of [0/45/-45/90] is a suitable one. However, in the present study as delamination is of interest and maximum interlaminar stresses are observed when maximum stiffness mismatch is achieved, i.e. maximum deviation in angles of adjacent plies is present, only the [0/90] architecture is considered. The methodology is nevertheless general and can accommodate any angle of orientation for fibres.
As for constitutive modelling of the RRL or cohesive layer damage initiation is stressbased but propagation is energy-based. In the simulation software program used (ABAQUS 2018) damage could be modelled using cohesive elements with continuum or direct tractionseparation formulations. This approach is similar to discrete crack propagation due to equality of J-integral with critical energy release rate. As in the current study the RRL is assumed very thin the variation of stress through the thickness is ignored and the layer can be modelled using the direct traction-separation cohesive model. The following equations are related to damage initiation and propagation in the context of cohesive element modelling.
In the direct traction-separation formulation damage initiation is dictated by the following hyper-surface in the space of non-dimensional tractions: Where 〈 〉 is equal to if > 0 and zero otherwise. Crack propagation is governed by modal energy release rates as follows: Where each numerator is obtained as the line integral of Eq. (15) geometrically depicted in   Since more often than not, failure occurs in a combination of modes the point on the hypersurface must be adjusted (see Fig. 22).  Table 3 shows the properties of the cohesive layer used in the analyses of this section.

Static analysis
The deformation (transformation) of the morphing blade subject to uniform lower surface vertical traction is shown in Fig. 23. Stress distribution is similar to the case of the stiffener, however, the level of stresses in the skins at the top and bottom are higher due to the fact that the hollow section requires higher stress at outer layers to compensate for the lack of material in the middle.  While it is obvious that there is deviation in the response from linear to nonlinear, it is observed that this deviation is greater for horizontal displacement component and occurs at a much lower level of displacement. This is similar to the case of stiffener discussed before. Fig. 25 depicts the resistance curve (force-displacement curve) corresponding to the vertical total force as a function of trailing edge vertical displacement. As compared to the stiffener, more pronounced deviation is observed when nonlinear geometrical effects are taken into account.
Study of static delamination is not considered for two reasons: (1) in order to have a realistic assessment of the failure in static mode nonlinear effects must be incorporated which requires implementation of an algorithm as Riks [54] or modified Riks [55]. These can only be implemented in an implicit code and cohesive elements are notoriously difficult to work with in standard finite elements as there are severe convergence issues subsequent to the initiation of the first crack and analyses usually get aborted, (2) even when a full picture is not required and the effect of geometric nonlinearity could be implemented intrinsically a quasistatic explicit analysis would be required. This analysis will take a long time to run and does not provide much information as it does not lead to the correct resistance curve to be used in a subsequent dynamic analysis [27]. It is furthermore not an integral part of a dynamic analysis thus lies outside of the scope of the present work.

Frequency analysis
Similar to the case of the stiffener a frequency analysis has been conducted to obtain the natural frequency of the first mode thus to determine the regime of response under pulse loading. As dynamic delamination is triggered based on the distribution of stresses which is, in its own right, dictated by the pattern of deformation in the blade, inclusion of higher modes (e.g. for impulsive case) gives rise to higher localised stresses where low stresses are normally encountered if only the first mode were to be considered (e.g. for quasi-static case). It is therefore helpful to acquire an insight into the regime of loading a priori so that we can expect different damage localisation modes. Fig. 26 shows the first two modes of the vibration for the morphing blade. The associated frequencies are f 1 =542 Hz (T 1 =1/f 1 =0.002sec) and f 2 =2516 Hz (T 2 =1/f 2 =0.0004sec). Similar discussion as for the stiffener in terms of nodes, antinodes, wave number, etc. holds here and in the interest of brevity this line of explanation is not repeated.

DELAMINATION IN MORPHING BLADE
Delamination in composite laminates has been the subject of extensive research both at structural and material levels with focus on inter-laminar shear and peeling stresses [16-18, 23,24,56]. Furthermore, dynamic crack propagation has been studied using 2D models (e.g. see [56]). However, the study of interfacial crack evolution in response to stress waves in a laminated structure is still a topic of active investigation form both modelling and experimental perspectives. Our aim here has been to develop a computational framework of predictive capability for the delamination under stress wave propagation. Dynamic fracture analyses are conducted using the cohesive layer model as the RRL.
Initially a single layer of cohesive was inserted to observe the possibility of static delamination. As mentioned before, it is possible, although time consuming and problematic, to attain results for static delamination. The result of static analysis conducted was reassuring as delamination in mode 2 was observed at the expected location. Fig. 27 shows the evolution of damage in the single cohesive layer inserted in the morphing blade. The concentration of shear stress at the edge gives rise to delamination as the damage parameter SDEG approaches 1 (SDEG=0 means intact and SDEG=1 means damaged or failed).
There will be relaxation of stress at the vicinity of this zone subsequent to failure and release of energy. Normal and shear interlaminar tractions alter based on the extent of damage and are calculated based on the damaged constitutive relations. 3.1. Analysis of dynamic delamination using cohesive elements As an initial remark, it must be mentioned that the actual spatial and temporal distributions of pulse loading are complicated and cannot be reduced to simple analytic functions. Even in the static case the pressure and shear tractions at a particular angle of attack can have the form represented schematically by Fig. 28.
The morphing blade has been decomposed layer-wise and at distances 10mm apart cohesive elements are introduced. This means roughly every 20 plies are fully tied and delamination is only allowed through these weak zones which is a rational choice as in any economic model one has to sacrifice a bit of accuracy for the sake of economy. To start, a pulse load of high intensity but short duration has been applied to the lower surface. The spatial distribution is uniform and the pulse temporal shape is rectangular of duration 0.1 msec. The amplitude of the shock has been 5MPa to allow damage to be experienced by the blade. Fig. 29 show the evolution of damage in the blade.
Different stages of deformation have been explained in figure captions. As it is observed for the impulse imparted by this shock loading the damage is excessive and severe. It is expected that the blade loses its entire functionality due to this scenario of loading.
Along the same line a pulse load of uniform vertical traction distribution of intensity 2MPa and duration 0.1 msec was applied to the lower surface. The temporal shape of the pulse was an isosceles triangle. Fig. 30 shows the excessive shear deformation of elements on the top of the adhesive layer adjacent to the hole as a geometric discontinuity under impulsive triangular loading. It is expected that delamination should start here, however, due to low amplitude of the loading no deamination occurs in this case.
Subsequently, a triangular pulse load of amplitude 2 MPa and duration of 2 msec was applied and delamination instigated at the expected position and the interlaminar crack propagated until total detachment of one of the arms and loss of functionality of the entire blade. Fig. 31 shows the different stages of deformation.
As mentioned before, due to maximum elastic mismatch between the adjacent layers ([0/90] fibre alignment architecture) delamination is more likely in this case when compared to other fibre arrangements. Fig. 36 shows the distribution of stress in the case of an impulsive load of amplitude 0.5 MPa and duration of 0.1 msec. The sharp change in stress from one ply to the next is obvious.
Finally, a quasi-static triangular load of amplitude 2MPa and duration of 20msec is applied. General failure pattern in this case is similar to the dynamic case, however, since the loading duration is higher and the load is more gradually applied there is chance for formation of multiple cracks as observed in Fig. 33.

Analysis of dynamic delamination using peridynamics
Here peridynamics is used to simulate fracture patterns and loads in the morphing composite blades. In particular through inclusion of intra-laminar failure this approach is suitable to test whether there would be any interaction between inter-and trans-laminar fractures in the blades modelled.
It may be expected that the use of XFEM in conjunction with the cohesive model may yield the same results, however, XFEM has not been used as anisotropic asymptotic crack tip functions are not implemented in the software ABAQUS and separate coding is required to implement these effects. An XFEM code developed at Imperial College London could be used in theory, however, two associated issues do not allow for this: 1) the XFEM code must be substantially modified to allow for cohesive interaction and accurate propagation of cracking, and 2)the XFEM code must allow for complicated geometries such as the morphing blade to be accommodated. These capabilities are not developed yet.
In a previous work, the bond-based peridynamic theory was extended to modeling fracture in anisotropic media by using spherical harmonics [49]. Here we used this method to simulate inter and translaminar fracture in the composite morphing blades. For the sake of comparison and since the capabilities allow the composite stiffener, has also been included as an alternative to the aluminium stiffener. For the aluminum stiffener, the isotropic formulation of the theory was used. The material properties used for the simulations were elastic modulus, density and fracture toughness as presented in Table 1. For the composite layers, elastic modulus and fracture toughness in fibre and matrix directions were used. To model interlaminar fracture, the average of the stiffness constants in fibre and matrix directions and a translaminar fracture toughness was used. Material properties are shown in Table 3 plus the intarply fracture energies assumed as G1=31.9 and G2=2.2 N/mm [53] which refer to G-values in the fibre and transverse directions, respectively and not to modes I (normal) and II (shear) fracture.
We used an explicit solver the details of which could be found in [49]. The stable time increment that we used was approximately 14ns. Simulations were performed on an Intel Xenon (3.5GHz) CPU using 16GM RAM. The simulation wall-time for 0.2ms duration of the simulations was approximately 3 hours and 45 minutes and linearly increased with loading duration. Fig. 34 shows the three models set up using the in-house developed peridynamic solver software for the different materials considered. In the composite blade (Fig. 35), the crack initiated between the composite layers (delamination) and continued as a delamination until the load bearing strength of the structure was significantly reduced. In this phase no trans-laminar crack was observed. However final catastrophic failure involved both inter and trans-laminar fracture. The location of crack initiation was different in the composite blade. The crack initiate near the bottom edge of the blade in the morphing composite blade and grew horizontally between the layers (Fig. 36). The strength of the blade were different with the aluminum blade being the strongest and the morphing the weakest thus the most vulnerable as well as the most flexible.

DISCUSSION AND CONCLUSIONS
The current study deals with presenting a computational modelling scheme for the study of morphing composite blades and their stiffeners which are used in aeronautical as well as green energy industries as an economically advantageous method for pitch control and maximisation of lift-to-drag ratio and delaying stall. This leads to optimising energy production in wind turbines and minimisation of fuel consumption for aircraft.
Aluminium stiffener and composite morphing blade are studied. Both extended finite element method (XFEM) and nonlocal continuum mechanics (peridynamics) have been used to study fracture in the stiffener. As for the composite morphing blade, cohesive elements are used to represent the RRL and delamination has been studied under shock, impulse, and pulse loads. The reason XFEM cannot be used to do so is the fact that asymptotic crack tip functions change for anisotropic materials [57] and have not been implemented in the commercial software used (ABAQUS 2018).
While both cohesive elements and peridynamics have been successful in capturing the pattern of delamination under uniform pressure and/or uniform lower surface vertical traction, under pulse loading only the results for cohesive elements are presented. This is due to the fact that peridynamic simulations will take a very long time and were merely conducted to conclude the similarity between the final results. Furthermore, finite element simulations using ABAQUS 2018 are much more commonplace in academia and industry and may be replicated by people who have access to the software while the peridynamic code is an inhouse developed code used by advanced users at Imperial College London.
The differences and similarities between delamination patterns for impulsive, dynamic, and quasi-static loadings are appreciated and in each case detailed analyses of delamination patterns are presented. It is anticipated that in the case of an impulse higher modes of vibration are triggered and the delamination occurs throughout the blade in several locations simultaneously or sequentially, however, for quasi-static loading the pattern of fracture is more as propagation from a single location till complete separation and subsequent cracks are mostly geometric i.e. due to large deformations and change of initial configuration. For dynamic pulse a mixture of both phenomena may be encountered depending on whether the pulse is closer to an impulse or to a quasi-static load.
While it is generally recognised that initiation of cracking is problematic the analyses are carried out to appreciate the degree and mode of damage in each case up to total failure, if this occurs. Parametric and sensitivity studies may be conducted to appreciate the degree and direction of influence of model parameters involved (such as lay-up architecture, type of material, or fibre orientation), nevertheless, these have not been included as they would elongate the paper needlessly and provide little further contribution to the main core objective.
As a remark, it must be mentioned at this point, that the actual spatial and temporal distributions of pulse loading are complicated and cannot be reduced to simple analytic functions. For instance, the static pressure load at a particular angle of attack can have the form represented schematically by Fig. 28. As for pulse loads a more realistic pattern of loading may require a fully coupled fluid-structure analysis, a subject which is outside the scope of the present work. Such studies usually involve a coupled arbitrary Largrangian-Eulerian formulation and a domain decomposition algorithm [13,[58][59][60].
In conclusion, it is shown the proposed scheme is capable of capturing all major features of dynamic delamination in composite morphing blades through an accurate, coherent, and consistent computational modelling scheme. This has a major impact in design applications where dynamic pulse and impact loads of all natures (accidental, extreme, service, etc.) are to be considered in design and therefore the framework may be used in design of lightweight morphing blades and wings.