Numerical and experimental investigation of out-of-plane fiber waviness on the mechanical properties of composite materials

The limited capability to predict material failure in composite materials and specifically in wavy composite layers has led to high margins of safety for the design of composite structures. Thus, the full lightweight potential of this class of materials is left unused. To understand the complex failure behavior of composite materials containing out-of-plane fiber waviness under compressive and tensile loading, a non-linear 2D material model was implemented in ABAQUS and validated with extensive experimental test data from compression and tensile tests. Each test was recorded by a stereo camera system for digital image correlation to resolve damage initiation and propagation in detail. This study has shown excellent agreement of numerical simulations with experimental data. In a virtual testing approach various parameters, i.e. amplitude, wavelength and laminate thickness, have been studied. It was found that the failure mode changed from delamination to kink shear band formation with increasing laminate thickness. The wavelength has shown minor influences compared to amplitude and laminate thickness.


Introduction
Fiber-reinforced composite materials allow for a significant reduction in weight due to the comparably low density (c.f. 4-5 times less than steel) and, in addition, fibers can be aligned in accordance with to the load paths. This possibility of alignment allows to place the fibers exactly at the position where they are needed to provide the component with its needed stiffness and strength. However, this can lead to a load path-optimized composite structure, not necessarily manufacturable in a robust and defect-free way.
The placement of the fibers or semi-finished textile products is still often carried out by hand, especially in the aviation industry. This allows diverse draping of the unidirectional (UD) layers, woven textiles or non-crimped fabrics (NCF) into the production tool. However, manufacturing effects such as fiber waviness, porosity, delamination and distortion cannot be completely avoided. The increased demand for composite parts for the aviation and automotive industries requires a transition to (partially) automated manufacturing processes. Those systems come with a higher deposition rate and ensure reproducible quality, but also imply production effects, e.g. fiber waviness [1,2]. This necessitates a well-founded understanding of those implicit effects on the mechanical properties of the manufactured structure.

Wave description and influence parameter
Fiber waviness is denoted as a wave-formed ply and/or fiber deviation from a straight alignment in a unidirectional laminate. Wavy plies can appear in arbitrary shapes and locations and can principally be classified into in-plane and out-of-plane waves, whereas Nelson et al. [3] stated, that both show similar strength degradations.
In general, the shape of fiber waviness is described by the ratio of amplitude to wavelength. Davidson and Waas [4] introduced 6 parameters to characterize the wave, indirectly including the laminate thickness. While amplitude A and wavelength L are most commonly used to describe the wave ( Fig.1), the maximum deviation θ max is considered to have the greatest impact on the mechanical properties.
Important influencing parameters such as maximum deviation of the fiber orientation from the original orientation, as well as the laminate thickness, have received little attention in previous investigations. The wave pattern is typically represented mathematically as sinusoidal waves [5,6]. El-Hajjar and Petersen [7] used a Gaussian function to capture the bell curve of wavy plies, which was found to better represent the wave geometry. Further, a more outward lying wavy layer leads to a greater decrease in strength, since the layers on the edge of the laminate are supported only one-sided and therefore fail earlier [8].
As the matrix properties and especially the thereby caused shear nonlinearity of the composite plays a dominant role in the failure behavior of wavy layers in composites, hot-wet conditioning (typically 70°C and 85% relative humidity in aviation industry) is also considered to substantially influence this behavior.
For a good prediction quality, the geometry of the fiber waviness must be reproduced as exactly as possible or in its stochastic deviation. Creighton et al. [9] and Sutcliffe et al. [10] have captured the fiber orientation and its scattering by means of digital image analysis. Nikishkov et al. [11] investigated automated FE mesh generation based on CT images.

Occurrence and assessment of fiber waviness
The origins of fiber waviness are manifold and a good overview is given in [1,2]. Some relevant sources are the draping step of the layers whether by hand or automated, compacting the layers in the autoclave or in the RTM tool [12,13], and gaps, overlaps, bridging and fiber steering in the AFP process [14,15].
Those effects cannot be completely avoided and therefore have to be tolerated and be considered as an inherent part of the structure. Fiber waviness can be considered as one of the most significant effects occurring in composite materials, due to the severe knockdown effect on the mechanical properties, such as stiffness, strength and fatigue and therefore, dramatically reduces the load carrying capacity of the material. If these unwanted irregularities are categorized as manufacturing features, i.e. effects, or defects, depends on the size, amount, and location of the effects in the component. On one side, the strength and stiffness reserve at the location of the feature and on the other side functional requirements, e.g. water tightness, have to be taken into account. The assessment of manufacturing effects further depends on the industry. In aircraft industry, the allowance limits for defects are very restricted, while in the automotive industry the requirement of short cycle times lead to a trade-off between robust processes and tolerated manufacturing imperfections.
To this point, there is still no generic acceptable approach to quantitatively support accept/reject/repair-decisions and make a consistent assessment of wavy layers in composites. If the feature is decided to be a defect, typically a deviation from design must be requested in aviation industry and an individual decision must be made on "use as is", repair or rejection. In some cases representative experiments on subcomponentlevel are performed on a statistical basis, however, this is both time consuming and cost intensive.
Therefore, it is necessary to strive for a fiber-, and especially a manufacturing-oriented design and construction of composite components. Towards this goal, design and production Fig. 1 Parameters influencing the mechanical properties, e.g. amplitude A, wavelength L, laminate thickness t, maximum fiber misalignment angle θ max , environmental conditions (room temperature or hot-wet), layup and the position of wavy layers within the laminate engineers aim to expand the permissible margin of safety by assessing the effect on stiffness and strength of those production effects, i.e. fiber waviness, porosity, delamination etc.. Additionally, they aim to reduce, or in the best case avoid, them on the process side, increasingly with the help of process simulations such as in [16][17][18].

Mechanical response of wavy composites
Even simple global loading cases, e.g. uniaxial extension, of wavy layers in a composite material locally induce a complex three-dimensional stress state consisting of a combination of interlaminar normal and shear stresses. Interlaminar tensile stresses normal to the lamina plane are induced as a consequence of transverse strains caused by the wave geometry [19]. Groundbreaking experimental and analytical studies on the influence of fiber waviness on mechanical properties were performed by Hsiao and Daniel [5,20] and Piggott [21].
It has been shown that compressive strength and stiffness are especially sensitive to the presence of fiber waviness. The initial misalignment can decrease the compressive strength by up to 60% [22].
Experimental studies under compression load for fiber waviness were carried out by [23][24][25][26], all reporting a significant reduction on compression properties. Lemanski and Sutcliffe [27] numerically investigated the influence of fiber waviness under compression load. It was shown that the compressive strength falls rapidly with the proportion of laminate width covered by the wavy region. The effect of length and position of the wavy region was found to have a smaller effect on the compressive strength. Zhu et al. [28] developed a threedimensional analytical method based on the classical lamination theory. Li et al. [29] have dealt with the determination of the reduced stiffness due to fiber waviness by proposing a micro-sphere based homogenization approach.
In a numerical and experimental study conducted by Vogler et al. [30] a strong interaction between compressive and shear stresses on compressive strength could be demonstrated. The general failure behavior of composites under compression and its resulting failure mechanisms, i.e. elastic microbuckling, plastic microbuckling, delaminations, fiber crushing and longitudinal cracking, have been studied in [4,[30][31][32][33][34][35][36][37][38], and similar mechanisms can be observed in composites containing fiber waviness. It has been shown that the nonlinear matrix properties dominate the damage initiation and the stiffness and strength decrease [22,39].
The formation of a shear kink-band is induced by local fiber buckling or matrix flow induced by initial or asdelivered fiber misalignments in the UD material, which continue to rotate under load. In the case of fiber waviness, this effect is further intensified due to the characteristics of the waviness. A very detailed study on experimental, analytical and numerical aspects on the formation of kink-bands can be found in [40,41]. Davidson and Waas [4] have shown that under compression load this is the critical failure mode also for laminates containing fiber waviness. They used a localglobal approach, where the local FE-based micromechanical analyses for fiber kinking are integrated into the global macroscopic FE model of fiber waviness, which exactly maps the geometry of the fiber waviness. The World-Wide Failure Exercise of Hinton et al. [42] have shown that the failure criteria show significant weaknesses, partly due to not considering fiber kinking under combined compression and shear stresses.
Under tensile loading, the fibers are pulled "straight", resulting in a stiffening effect and a failure of the interface between fiber and matrix due to local interlaminar normal and shear stresses and finally to delamination [43,44].
For shear loading in the laminate plane of a laminate with fiber waviness, few experimental and simulative investigations are known. In the few available literature references [15,45] an influence of up to 12% on the reduction of the shear strength is determined. Similar to shear stresses in laminates without fiber waviness [46], a strong influence of the nonlinear shear behavior on the failure process is assumed for laminates containing fiber waviness.
In recent publications the simulations are validated on strain field measurement to record damage initiation and progression [4,47]. Previous studies have mainly focused on UD laminates, all oriented in the same direction, and have concentrated on the effect of fiber waviness on the stiffness. Quasiisotropic (QI) laminates containing wavy layers have been rarely investigated so far. Mukhopadhyay et al. published work on compressive [48], tensile [49] and fatigue behavior [50] of QI laminates containing embedded wrinkles using 3D finite element modelling approaches reporting a more complex failure behavior due to multiple active damage mechanisms and their interaction.
The study of Wilhelmsson et al. [51] has shown a considerable influence of embedded out-of-plane waviness on the bending during compression testing.
This study aims to investigate the mechanical behavior of out-of-plane fiber waviness of both unidirectional and quasiisotropic laminates by numerically simulating damage initiation and propagation, using a custom material model, implemented in ABAQUS/Explicit [52]. To determine material parameter and validate our simulations we conducted uniaxial compression and tensile tests supported by optical strain field measurements.

Constitutive material behavior
The constitutive material behavior of a single UD ply can be modelled by an undamaged non-linear elasto-plastic hardening regime and an elasto-damage softening regime shown in Fig. 2. In the hardening regime the stress increases with ascending strain. The degree of non-linearity depends on the angle between loading direction and orientation of the layer. The failure criterion according to Puck [53] is used to determine the damage initiation. The degradation of the continuum is described in the softening regime and controlled by energybased linear stiffness degradation. Modelling of interface delamination, e.g. cohesive zone modelling (CZM) or virtual crack closure technique (VCCT), is not part of the model. For laminates with all plies in 0°-direction the use of a 2D model in combination with a plain strain assumption has been validated in [22]. Testing the applicability of this approach for QI laminates is still an open task and part of this current study.
Non-linearity in the elasto-plastic hardening regime A plasticity model developed by Sun and Chen [54] is implemented to describe the non-linear constitutive behavior of a IM7 unidirectional carbon fiber reinforcement embedded in a 8552 thermosetting epoxy matrix. The plasticity model (plastic potential) is based on a quadratic yield function f(σ ij ), defined in terms of the stress components σ ij in the principal material directions.
The plastic potential is expressed by the general yield function where the material constant a 66 describes the anisotropy in the plasticity.
The incremental strain dε ij can be decomposed into an elastic strain component dε e ij and a plastic strain component dε p The in-plane incremental plastic strain components dε p ij are defined in terms of the in-plane yield function f and the plastic multiplier dλ via the associated flow rule. Assuming that the longitudinal plastic strain increment dε p 11 can be neglected since most unidirectional composite materials behave almost linear elastic in longitudinal tension and compression loading [54] the plastic strain component can be described as where γ 13 denotes the engineering shear strain. Similarly, the work per unit volume dW performed during elongation dε ij can be also decomposed into a reversible elastic work dW e and an irreversible plastic work dW p . The incremental plastic work dW p is defined as By introducing the effective stress σ ¼ ffiffiffiffiffi ffi 3f p for the case of a plane stress state, the plastic work increment dW p can be written as Rewriting in terms of the effective plastic strain increment dε p leads to The incremental plastic stress-strain relation only depends on the values of the plasticity coefficient a 66 and the plastic multiplier dλ which is calculated as with the plastic modulus H p H p is a derivative of the effective stress with respect to the effective plastic strain and using the current effective plastic strain.
Fiber composites are assumed to have no clear yield locus. For this reason, a power law approach can be used to describe the non-linear constitutive behavior in order to fit a master curve of equivalent stress over equivalent strain. The compression and tension regimes must be considered independently of each other.
For a complete description of the non-linear material behavior, the parameter a 66 must be obtained from off-axis tests (OAT).
In the plane stress state, the prevailing stresses can be indicated as follows with the associated strains Assuming a purely elastic material behavior, stress and strain can be related by means of Hook's law (linear elastic material law) by introducing the material stiffness matrix, respectively compliance matrix. In Voigt notation in combination with Einstein's summation convention follows.
For the plane stress state, the stiffness matrix C 0 ij can be defined as The corresponding compliance matrix S 0 ij is written as : ð14Þ Continuum damage model

Damage variables
In the case of damage progression in the lamina, the constitutive relationship resulting from the progressive damage is considered in the form of a reduction of the stiffness matrix. In this case, a pure elasto-damage is assumed. A damage variable d was introduced in the form of a thermodynamic state variable to determine the damage progression. In the general anisotropic case, the damage of fiber composites is described in the form of a 4th order damage tensor D. In addition, a 4th order damage influence tensor M(D) is introduced to characterize the damage condition. Applied to the concept of the equivalence of stresses, the effective stresses in the damaged continuum can thus be written as wherein I represents the 4th order unit tensor. In the case of an uncoupled damage and a unilateral damage propagation in the direction that tears equivalent to the fiber or matrix direction, the damage tensor D can be written for the pane stress state wherein 1 and 3 denote the fiber and matrix direction, respectively. This allows the damage influence tensor M(D) to be specified.
In the case of material damage in the continuum mechanics, the damage variable d is introduced as a thermodynamic state variable. The evolution of the damaged stiffness matrix by internal state variables d are schematically shown in Fig. 3 for transverse direction. The free energy potential f d can be defined as a function of d. For the damage process free of thermal influence, the following should therefore apply The strain equivalence approach can be used to define It finally follows for the damaged stiffness matrix C d (D) ð23Þ and the elasto-damaged material law The damaged compliance matrix S d (D) results equivalent to The dissipation inequality for the damaged material behavior results in two conditions. According to Coleman and Gurtin [55], it applies as for the undamaged material Thus, the following condition can be derived from the dissipation inequality, which must be fulfilled for the inequality to apply The term ∂ f d ∂D corresponds to the energy release as a result of progressive damage to the material. This term is called energy release rate Y and an energetically conjugated quantity of the damage variable.
Y Ḋ describes the energy dissipation rate in the material, where Y is the thermodynamic force and Ḋ is the thermodynamic flow.

Damage initiation
For the damage initiation the well-established Puck criterion [53] is used. The individual failure modes defined by Puck determine the progress of the damage.
Fiber failure under tension σ 1 ≥ 0 Fiber failure under compression Inter-fiber failure under tension Inter-fiber failure under compression Inter-fiber failure under compression σ 3 < 0 (Mode C) The linear approach for the constitutive relationship of the material in the softening regime results in where K characterizes the descent in the softening regime. The damage evolution law is based on the concept of effective stresses in the context of strain equivalence as shown in Fig. 4. The approach applies to uniaxial stress, where f Ei is the stress exposure related to a stress state in the undamaged continua.
This definition is used to determine the damage progression as a function of the effective stresses.

Energy consideration and element size regularization
In continuum damage approaches the damage is represented non-locally, i.e. related to a unit volume (continuum) homogenized in the constitutive relationship. Due to the brittle fracture behavior of fiber composites, the formation of macrocracks due to diffuse micro-damages takes place in a relatively narrow specific area known as the crack process zone (CPZ). In order to enable a relationship between the non-local definition of the continuum damage approach and the local definition of crack formation, the following relationship can be established. The product of the energy density per unit volume with the specific width of the crack process zone corresponds to the fracture energy resulting from the formation of the stress-free crack. In other words, the dissipated energy for the complete deterioration of the material in terms of a unit volume is equal to the dissipated energy for crack development in terms of a unit area.
Since the quantities in this equation represent all intrinsic material parameters, the numerical FE consideration would have to be discretized in such a way that the element size corresponds to the intrinsic width of the crack process zone. In the consideration of damage processes of fiber composite materials this procedure is not useful. Therefore, a new meshrelated width parameter l * is introduced, which is no longer bound to the physical width of the crack process zone, but represents a regularized numerically defined quantity. By adjusting the area under the stress-strain curve, a dissipated energy in the failed element is obtained as a function of the ratio of physical crack zone width and numerical characteristic element length.
An element size regularization is thus carried out by adjusting the descent in the softening area.
Both the strength and stiffness are independent material parameters. The fracture energy serves as an intrinsic input parameter. Mode I and II fracture toughness can be empirically determined by material tests using a DCB and ENF test, respectively [56].
Two configurations of sinusoidal waves were realized using a one-side female metal plate tooling in which the defined sinusoidal wave configurations where milled in with a wavelength L 1 = 27.9 mm and amplitude A 1 = 1.19 mm (wave 1) and L 2 = 14.5 mm and A 2 = 0.58 mm (wave 2). Laminate quality was verified by measuring the fiber volume fraction (61.3%, 0.95% STD) using wet chemical fiber extraction according to EN 2564 for 12 specimens.

Tensile and compression testing
Tensile and compression tests are carried out in the spirit of ASTM D3039 and ASTM D6641. Material testing was conducted on universal test machines (Zwick, Ulm, Germany) with a maximum load of 150 kN. The tests were displacement controlled with a rate of 2 mm/min and continued until failure occurred. An overview of the conducted tests is given in Table 1.

Digital image correlation
Digital image correlation (DIC) is a powerful tool to measure the full field strain distributions during mechanical testing of materials. DIC provided information on mechanisms and a full 2D strain field prior to failure, i.e. the location of strain concentrations. However, capturing damage initiation is limited by the (typically lower) frame-rate as well as the detectable size of an event. We used a stereo camera system with two 12 Megapixel cameras (GOM Aramis™ 3D 12 M) equipped with 100 mm lenses. The DIC measurement system was triggered by the test software via analog inputs, ensuring synchronized recording of images with the force and displacement signal. Images were obtained at 10 Hz. A homogeneously distributed speckle pattern was applied on a matt white grounding using the airbrush system Minijet 4400 B RP (SATA, Kornwestheim, Germany). The test-set up is shown in Fig. 5. To analyze the images we used GOM Correlate Professional (2016). The facet size was set to (14 × 14) pixels with an overlap to get a point distance of 10 pixels.

FE model
The finite element model of both the 15°off-axis tests and wavy geometries are modeled with reduced integration plain strain elements (CPE4R) including hourglass control to avoid zero energy modes. For parameter studies a python script was written to automatically build up FE models with varying parameters, e.g. amplitude, wavelength and laminate thickness. The material model described in Section 2 was implemented as a material usersubroutine VUMAT in ABAQUS/Explicit [52]. A representative FE model of a wavy specimen including boundary conditions and loading is shown in Fig. 6. For better visualization, a coarse mesh is shown. Prior to the numerical simulations, a mesh size study was carried out and 35-40 elements in thickness direction were found to be sufficient. Simulation results of the two wave configurations (wave 1 and 2) were compared with experimental test data. A parameter study was carried out for varying amplitudes (0.1A-2A), laminate thicknesses (0.5A-6A), where A denotes the baseline amplitude of 1.19 mm, and wavelengths (13.45-41.35 mm). For modelling quasi-isotropic layups with the 2D approach the material properties were transformed for each orientation according to classical laminate theory and applied to the corresponding layers.
The material properties of IM-7-8552 carbon fiber reinforced epoxy matrix prepreg used for the simulations are given in Table 2. Most simulations were carried out before experimental tests, therefore material parameter from literature were used. These values are consistent with our tested values (c.f. Table 2 and Table 3).

Results and discussion
Mechanical behavior of planar reference and wavy specimen The experimental data for planar reference and wavy specimen under compression (Figs. 7 and 8) and tensile loading (Fig. 9) show a non-linear stress strain response for wavy specimen and a considerable knockdown. Young's modulus in fiber direction and the peak-stress before failure for all tests are shown in Table 3. The compressive strength is reduced by approximately 75% for the more pronounced wave configuration (wave 1) and 50% for the less pronounced wave (wave 2). The latter also shows a substantial reduction of the strainto-failure by approximately 45%. Those results are in good accordance with previously reported findings [22]. Although the compression properties were assumed to be more affected by fiber waviness compared to tensile loading [43,44], results from mechanical tests on wave 1 show a drop of 70% in tensile strength for first ply failure, where mode I delamination occur due to positive stresses in thickness direction σ z (Fig. 13c). We found a dramatic drop in compressive strength for both wave configurations with quasi-isotropic laminate compared to the reference specimen. The compressive strength dropped by approximately 65% (wave 1) and 32% (wave 2). Interestingly, we could only observe a reduction in stiffness for the more pronounced wave (wave 1) but not the other one (wave 2). In compression loaded laminates, containing fiber waviness, the mechanical properties of the 90°-layers are not reduced by the wave. The highly affected 0°-layers  account only for 25% of the total laminate. Therefore, the presence of waviness has a reduced influence on QI laminates compared to UD laminate. This results suggest, that the drop in stiffness properties for the more pronounced wave 1 is mainly caused by the geometrical deviation of the wave due to a higher amplitude-to-thickness ratio (A/t). The resulting Young's modulus and strengths for quasi-isotropic laminates are summarized in Table 3.
The applied unidirectional load leads to a complex 2D strain state at the observed cross-section due to the presence of the wave. In Fig. 10 the shear strain ε xz distributions of wavy specimens are shown right before and after final damage within two subsequent frames (10 Hz) for a) wave 1 and b) wave 2 under compression and c) wave 1 under tensile load. Regions of high strain coincide well with the occurrence of macroscopic material failure. This is in accordance with [20], who determined that the associated shear stress τ 13 is the most significant stress component for damage initiation in wavy composites under axial compression.
While compression loads increase the amplitude-towavelength ratio, global applied tensile loads straighten the fiber waviness, leading to an increasing stiffness during loading [57]. Both of which require a geometrically nonlinear calculation. The substantial influence of both the amplitude A and thickness t in light of the resulting bending stresses are shown in Eqn. 40. A simplified analytical consideration using the Euler-Bernoulli beam theory (neutral axis coincides with the beam centroid line, bending moment M B = FA = σ x b t A, where A is the amplitude, the moment of inertia for rectangular cross section J y ¼ b t 3 12 and y ¼ t 2 ) leads to Thus, the resulting bending stresses, with a maximum in the longitudinal center of the wave, increases with increasing amplitude and decreases with increasing thickness. To incorporate the geometric non-linearity the curved beam theory can Fig. 9 Tensile test of UD-laminate -Comparison planar vs. wavy specimen Fig. 8 Compression test of QI-laminate -Comparison planar vs. wavy specimen Fig. 7 Compression test of UD-laminate -Comparison planar vs. wavy specimen be utilized, however a comprehensive analytical study is outside the scope of this contribution. The resulting shear stress distribution from a linear-elastic but geometrical non-linear FE simulation evaluated along a centroid curve for varying laminate thicknesses, shown in Fig. 11, underlines the thickness-dependent behavior. Shear stresses are highest at the location of the biggest incline/ decline and generally increase with thickness. However, when a thickness of 4A is reached, the further increase of the resulting shear stress is almost zero. An additional linear elastic FE simulation of varying amplitudes and constant wavelength L = 27.9 mm and laminate thickness t = 7.14 mm shows the superposition of bending stresses with the global compressive or tensile load in Fig. 12 by plotting the stresses in x-direction at the peak point of the wavy specimen. This results in different stress states (tension or compression) in the center of the wavy region, depending on the thickness. The compression loaded wave with an amplitude of A = 0.595 mm shows approximately zero stress σ x at the peak point due to the superposition of global loading and resulting bending stresses. This mechanism in form of variations of the longitudinal strain component ε x is also schematically shown in Fig. 13 and can be clearly observed in DIC measurements. There is a transition from tension to compression of ε x across the thickness when global compression and tension is applied within the tested wave with the higher amplitude (wave 1, c.f. Fig. 13a). Wave 2 (less pronounced waviness) shows a gradient in the strain component but no transition from tensile to compression strains. The effect of superposing bending stresses is amplified during compression testing due to an increasing amplitude-to-wavelength ratio (A/L).
In addition to the above mentioned shear strains that are influencing the failure behavior considerably, resulting tensile strains ε z , respectively stresses in thickness direction σ z potentially cause the individual layers to lift off, causing delamination, as shown in Fig. 13a and c). This leads to the conclusion, that the resulting stresses superpose with shear stress components and also considerably influence the failure behavior. In global tensile loading (Fig. 13c), the resulting through thickness strain component ε z in wave 1 shows a maximum in tension at the center of the wave, leading to delamination along the whole wave length, potentially initiated at the wave center. Wave 2 shows a maximum tensile strain component ε z at the turning point of the sinusoidal wave under compression loading. This leads to one-sided delamination initiated at the area of maximum inclination of the sine wave. Global compression load of wave 2 and the absence of notable bending moments leads to fiber kinking of the specimen schematically depicted in Fig. 13b). The positive strain ε z in thickness direction leads to less supported fibers and therefore promote the appearance of kink bands. These results suggest that the failure of wave 1 is determined by the bending moments resulting from the geometry (mainly amplitude and thickness). In contrast, for wave 2, the contribution of bending stresses in superposition with the global load is less thus we observe material failure. The failure behavior is more determined by the material properties for lower amplitude-to-thickness ratios.

Model validation
Off-axis tension tests are carried out to determine material parameter needed to describe the nonlinear behavior in combined transverse tension and in-plane shear. According to Körber [58], 15°off-axis specimens are most suitable for the determination of shear properties. In off-axis tests with straight end tabs, it was found that the extension-shear coupling leads to significant stress concentrations near the clamping regions and to a non-uniform strain field. To avoid this, oblique end-tabs were used as proposed by Sun and Chung [59]. The stress-strain response from both simulation and test data for 15°off-axis tensile tests is shown in Fig. 14.
The numerical model is able to match the nonlinear material behavior and the damage initiation properly. However, the resulting failure strain/stress was typically too low, thus our model serves as a conservative estimate. Remote stress-strain curves stemming from simulation and experimental results are shown for wave 1 and wave 2 under compressive loading in Figs. 15 and 16, respectively. The shear strain distribution for wave 1 shows excellent agreement of numerical and experimental results. The global stress-strain relationship for wave 2, including the onset of failure, location and failure pattern shows excellent agreement of simulation and test data.
The typical failure mechanism of a narrow kink-band formation due to matrix yield under compressive loading are indicated at locations of maximum inclinations. The nonlinear stress-strain response is increased at higher remote strains. In Fig. 15 a slightly conservative damage initiation at approximately 50 MPa below test results is probably due to the use of stress-based failure criteria. When the stressstrain curve becomes strongly nonlinear and reaching a plateau, the exact point of damage initiation is more difficult to find. A slight increase in stress leads to a considerable increase in elongation, thus strain-based criteria might be advantageous. The numerical model result is compared with experimental results, determined by digital image correlation (DIC), by comparing the shear strain distribution along the center line of the wave in Fig. 17. The axial strain ε x from our simulations and the DIC measurements evaluated along a centered path vertical to the loading direction shows good agreement of numerical and experimental data for both wave configurations with UD laminate in Figs. 18 and 19.
To evaluate the applicability and limitations of the implemented 2D continuum damage model for more complex materials, it is applied for wave 2 under compressive loading for a quasi-isotropic laminate (0/45/−45/90) 5s . The global stressstrain response, shown in Fig. 20, can only be reproduced insufficiently. The resulting complex 3D stress state leads to a complex failure behavior of QI laminates that cannot be represented by the 2D model, due to neglecting the stress components in y-direction. However, the resulting failure mode as well as the correct location for damage onset were detected correctly at the center of the wave. For the simulation of QI laminates, the model can merely be used as an engineering estimation of the material properties and a qualitative assessment of failure. Therefore, for analyzing QI laminates containing fiber waviness an extension to a threedimensional model which incorporates nonlinearity parameter for various ply orientations, fiber rotation and cohesive layers for modeling delamination is suggested.

Virtual testing
The model has been successfully validated for UD materials. To reduce the effort of physical testing the model is applied for studying the influence of varying parameters, e.g. amplitudes A, thicknesses t and wavelengths L, in the sense of a virtual test approach. The considerable influence of laminate thickness on the failure behavior has already been described in Section 4.1. The virtual testing   Fig. 21, illustrate the shift from failure dominated by bending stresses to the formation of kink bands with increasing laminate thickness. This confirms our statement, that for lower A/t ratios the testing regime represents a test at material level, whereas the failure behavior above a certain threshold for A/t ratios is strongly influenced by introduced bending moments, thus representing a test on structural level. This results in a pronounced non-linearity in the stress-strain curve of thin laminates (e.g. t = 0.5A and t = 1A) due to the influence of the geometry (geometrical non-linearity). Figure 22 shows the global stress-strain response of compression loaded waviness with varying amplitudes and constant thickness (t = 7.14 mm) and wavelength (L = 27.9 mm). Even minor fiber deviations associated with the amplitude lead to a drastic reduction of the mechanical properties, which is known from transformation relations for engineering constants as a function of different fiber orientations [60]. According to that, simulation results show the relative reduction of the stiffness to be smaller at higher amplitudes.
The results from FE simulations of waves with varying wavelength and constant amplitude (A = 1.19 mm) and thickness (t = 7.14 mm), shown in Fig. 23. The relative drop in compressive strength from L = 41.35 mm to L = 13.45 mm is~46%. However, compared to a planar reference specimen those specimens show a reduction of 76% (L = 41.35 mm) and 87% (L = 13.45 mm, thus only a 11% absolute difference. This confirms to Lemanski and Sutcliffe [27], who stated that the effect of defect length is relatively small in this configuration.

Conclusions and outlook
We successfully tested and simulated the mechanical behavior and damage initiation for two wave configurations with UD and QI laminates under compressive and tensile loading. The understanding of intrinsic material behavior of fiber reinforced composite materials on microscopic level up to macroscopic or structural level is of crucial importance for the development of suitable material laws for numerical modeling and for a deeper understanding of deformation and damage mechanisms. On top of that, understanding the material behavior of wavy composites is of vital importance for further cost and weight savings. The implemented nonlinear material model was able to capture the material behavior including shear nonlinearities, failure initiation and propagation in unidirectional  Fig. 23 Stress-strain relationship for wavy specimen (A = 1.19 mm, t = 7.14 mm) with UD laminate and varying wavelength using implemented non-linear VUMAT (Puck) Fig. 22 Stress-strain relationship for wavy specimen (t = 7.14 mm, L = 27.9 mm) with UD laminate and varying amplitude using implemented non-linear VUMAT (Puck) laminates reasonably accurate. Overall, our results suggest that, depending on the laminate configuration (amplitude, wavelength, thickness) inter-laminar shear failure is the dominant failure mechanism followed by mode I delamination and layer-wise buckling. The nonlinear shear behavior under compressive loading increases with the severity of present fiber waviness and, therefore, plays a dominant role for the initiation and propagation of the failure mechanisms. For wave configurations with a lower amplitude-to-thickness ratio fiber kinking is the dominant failure mechanism under compression loading. The observed formation of a shear kink-band under compression and shear stress (Fig. 10b) is caused by local fiber buckling, induced by misaligned fibers in the UD materials that continue to rotate under load. In this work, the significant influence of laminate thickness was shown, which has been rarely considered in both industrial practice and academic research so far. The failure behavior for lower amplitude-to-thickness ratios is more determined by the material properties, whereas higher amplitude-tothickness ratios lead to a higher influence of the geometry due to bending. For this reason, we suggest to consider the amplitude-to-thickness ratio (A/t) in addition to the commonly used ratio of amplitude-to-wavelength (A/L) for further assessment strategies of fiber waviness. The validity of the model applies for quasi-static loaded unidirectional laminates only. For multidirectional laminates an extension from 2D to 3D is necessary to incorporate the complex 3D stress state.