Modeling Shear Behavior of Woven Fabric Thermoplastic Composites for Crash Simulations

Woven fabric thermoplastic composites possess high specific strength and stiffness along with thermoformability. To utilize the full potential of these materials to achieve better crash-safe designs in automotive structural parts, the measurement of non-linear shear behavior and its material modeling for FEM simulations is required. The standard testing method was used to measure the pure shear behavior of woven fabric composites. These results were compared with the shear behavior of material in the presence of normal stresses along the fiber direction. Tensile and compression cyclic testing of ± 45° laminate were carried out to measure the stiffness degradation and hardening of the material in the presence of tensile normal and compression normal stress. A methodology is proposed for taking into account the differences in shear behavior under different loading directions in an FEM simulation. Based on the experimental evidence, improvements in the mathematical description of plasticity and damage in continuum damage mechanics models are proposed. The model was implemented as a user-defined material subroutine (VUMAT) for Abaqus. The experimental results from coupon tests were used to verify the results of a single element simulation. Finally, a three-point bending test was used to validate the predictions of the user material model.


Introduction
Woven fabric composites are increasingly being used in the automotive sector due to their high specific strength, specific stiffness, and specific energy absorption capabilities. The woven fabric thermoplastic composites, so-called Organo sheets, can be thermoformed, which makes them particularly feasible for series production with reduced production costs. To use these materials in automotive applications, crashworthiness requirements must be met.
The lightweight potential of these materials is achieved by aligning the fibers in the direction of the load; this, however, leads to catastrophic brittle fracture at low strains. Sometimes, it might be desirable to achieve ductile or high failure strain behavior which can be achieved by angle-ply laminates. Furthermore, a weight-efficient design in which fibers are aligned along one particular direction might not be in the direction of crash forces. Therefore, the ability to predict the crash behavior of angle-ply laminates through FEM simulations, then, is just as important as it is along the fiber direction.
Angle-ply laminates demonstrate non-linear behavior and high strains to failure. This effect is more pronounced in the tensile/compression loading of 45° off-axis laminates. The composite material response is linear elastic for 0°, and non-linearity increases gradually with the angle of loading. The source of metal-like pseudo-ductility in angle-ply laminates is the non-linearity of the shear behavior of composites. The shear stress-to-applied stress ratio is highest for the loading angle of 45° [1][2][3]. Therefore, testing and modeling of the shear behavior of composite material is essential in the product development process.
Woven fabric composites exhibit higher stiffness and strength under shear loading in comparison to unidirectional composites, provided that the fiber-matrix type and fiber volume ratio remain the same [4]. Measuring the shear properties of composite materials has been standardized through the use of various test methods [5][6][7][8][9][10].
Several orthotropic material models accounting for non-linearity in shear have been implemented in commercial FE code Ls-Dyna [11]. The *MAT_58 uses a bilinear curve defined by two pairs of shear stress-strain values to describe the non-linearity in the shear behavior of laminated fabric composites. A similar approach was adopted by Pinho et al. [12] by inputting the shear curve of laminated fiber composites into an FE model using a polynomial equation. The coefficients of the polynomial were calibrated to represent the value of shear stress as a function of shear strain. In this approach, the value of the shear modulus remains constant as represented by one of the coefficients in the polynomial. However, such approaches do not account for plasticity and damage in contrast to experimental data in which the degradation of shear stiffness due to damage is evident.
Cousigné et al. [13] used the Ramgerg-Osgood model for prediction of the non-linear behavior of material and a unidimensional plasticity model in a separate subroutine to account for permanent deformation. The model is used to calculate plasticity when a difference occurs between the initial and instantaneous stiffness. However, none of the four proposed damage formulations (linear, non-linear, constant stress level, step-based damage) are in agreement with the measured test data of organo sheets. Many authors [1,[14][15][16] have shown that the effect of fiber rotation is also essential for predicting the results. However, at large strains, the fiber rotation calculation for woven fabric is much more complex due to the curvature and locking effect in warp and weft, for which a micromechanical model must be used. This, in turn, will be computationally very expensive for crash simulations. Van Paepegem, et al. [17,18] carried out an experimental program to measure the shear properties of glass fiber-reinforced epoxy composites with tensile tests on ± 45° and 10° off-axis composite laminates. The results were used for establishing a phenomenological model of shear damage and permanent shear strains. Validation simulations of the three-point bending of rectangular ± 45° laminate have shown that measuring shear parameters only in the presence of tensile normal strains was not enough. Instead, shear parameters must also be determined in the presence of normal compression strains.
Ladeveze and LeDantec [19] were the first to use the Continuum Damage Mechanics (CDM) approach to model the non-linear shear behavior of uni-directional composite ply material, coupling it with plasticity. This formulation takes into account the damage by use of damage variables and their associated thermodynamic forces, effective stress levels, and elastic and inelastic strains. A damage variable is associated with the formation and evolution of microcracks and cavities. These defects are irreversible and can cause stiffness degradation. PAM-CRASH composite model 131 offers various damage evolution functions, including linear form as well as exponential form, between damage and the damage-driving force. However, only one set of parameters is used to determine the shear curve, regardless of tensile or compression loading. Johnson [20] generalized the concept of Ladeveze's [19] for the modeling of a woven composite. The classical plasticity with an elastic domain function and hardening law was applied to effective stresses in damaged material. However, the damage evolution function is not consistent with measured experimental data. Furthermore, the effect of tensile or compression loading in combination with shear loading at large strains was ignored.
The focus of this study was on testing and modeling of the shear behavior of woven fabric thermoplastic composites. Therefore, the state-of-the-art CDM model proposed by Johnson [20] served as the basis for this investigation. Major gaps were identified in state of the art material model and experimental data. State of the art material models ignored the change of shear behavior in presence of tensile and compression normal stress. And shear damage evolution function did not comply with experimental data. This paper addressed these problems and improved material model is proposed. Better simulations results are obtained by using improved material model.
In Sect. 2, a brief overview of the state-of-the-art CDM model proposed by Johnson [20,21] is given. In Sect. 3, the investigated woven fabric composite material is described, and in Sect. 4 an explanation is provided as to how established testing techniques are used to measure the shear stress-strain curve of the material. Experimental measurements of (1) pure shear behavior and shear behavior in the presence of (2) tensile and (3) compression normal stresses were carried out, using three different test methods and improvements which are proposed from the perspective of experimental evidence. In Sect. 5, improvements to the CDM model based on the experimental evidence are proposed. The extended CDM model is then implemented as a user defined material model (VUMAT) in Abaqus, and the working ability of VUMAT is verified on a single element simulation. Finally, in Sect. 6, the validation of simulation results on the three-point bending of ± 45° laminate is included.

Continuum Damage Mechanics
From the perspective of transient dynamic simulations, finite element-based codes with explicit formulations are used. Explicit formulation has the advantage of treating instabilities and material softening, especially with problems involving contact. For a reasonable computation time, multi-layered shell elements are used to represent plies [22]. The material modeling is possible at the multiscale or mesoscale levels. In multiscale modeling, a constitutive model is defined at the microscale (fiber and fiber-matrix interface level) and the resulting stress/strain fields are transformed to the meso-level. However, this form of modeling requires huge computational resources and is undesirable for automotive crash simulations. Therefore, mesoscale modeling is appropriate where lamina is considered as homogenized material [23]. CDM is a mesoscale model and considers the orthotropic damage elastic response of the material at the homogenized level. In CDM [19][20][21]24] the 1 3 intra-laminar material behavior of composite ply is written in material coordinates, keeping in view the plane stress state as: in which d 1 and d 2 are the damage along the fiber direction due to microcracks in fibers and d 12 is the shear damage due to matrix microcracking. The damage parameters take values between 0 ≤ d i ≤ 1 . The damaged-material strain energy in plane-stress state for homogeneous ply is written according to [19,20]: The so-called thermodynamic forces can be derived by taking the differential of strain energy potential with respect to damage variables. Considering a shear-dominated loading, in which material response along the fiber direction remains elastic, the shear stress can be written as [21]:

Plasticity
The stress-strain curve is linear below a certain value called yield stress. The stress in the material is bound by yield stress and constrained to the elastic domain. Before the yield point is reached, the strain is within the elastic limit and hence the strain is elastic. After the yield point, hardening takes place and the plastic strain increases. It is important to know that even after reaching the initial yield point, elastic strain can still increase. The elastic strain, then, can be thought of as the amount of strain that can be reversed during unloading. All the remaining strain after unloading is plastic strain. Therefore, the total strain is the sum of both the elastic and plastic strains, which can be written as [21]: Equation (3) can be rewritten to define the elastic relationship for effective shear stress according to [21]: Since the absolute value of the stress in material cannot be greater than yield stress, it can be mathematically described as [21]: in which F is called yield function. For every strain increment, the stress value is calculated and it must fulfill the yield condition. ∼ y is the yield stress due to hardening or the resistance of the material to deformation. In the plastic region, the resistance of material against the deformation is less than in the elastic region. The hardening starts after reaching the yield point, and post-yield behavior is seen as both increasing in elastic as well as plastic strain. When hardening takes place, the yield stress changes. The type of hardening after the yield point depends on the type of material. Ramberg-Osgood hardening curves are best suited for modeling the shear plasticity of woven fabric composites [21]: in which ∼ y0 is the initial yield stress and K and p are material parameters which can be identified from Fig.13.. | | | is higher than the yield value, which is not allowed. In such a case, the amount of plastic flow Δ is calculated such that F = 0 . This means that plastic flow remains zero within the elastic domain and nonzero on the boundaries which can be described mathematically as [25]: Now we define the term Δ mathematically in the form of the flow rule [21,25]: The sign function return -1 if the value of ∼ 12 is negative, and + 1 if the value of ∼ 12 is positive. In this equation, the term Δ is the change in the absolute value of plastic strain during a time step. Although it appears as a time derivative in equation, in an explicit simulation it is analogous to the plastic strain increment in a load step.
Plastic flow is calculated from the consistency condition ΔḞ = 0 [25]. In the case of Ramberg-Osgood hardening, equation Δ must be solved with the Newton-Raphson method to determine the increment in plastic strain [25].

Damage
It is understood that material cannot heal itself after damage has occurred. Therefore, after unloading, the damage cannot be reversed. The state of damage remains unchanged within a domain, defined by the damage activation function F 12 [21,24]: in which r 12 is the damage threshold, defining the range of the elastic domain. The initial value of the damage threshold is set equal to 1. 12 defines the damage initiation criterion in the form [21]: in which S is the damage onset stress value, measured experimentally. Since damage is a non-decreasing value, it must be monotonically increasing when damage takes place ( ḋ 12 ≥ 0 ). The damage threshold is not an independent variable, as it relates to the damage variable. It is defined by [21,24]: Therefore, the value of r 12 at time t is equal to the maximum value of 12 during all the history of time . This is important for implementation of the irreversible damage behavior of the material. What remains now is the definition of the shear damage variable, which was proposed by Johnson [20] as: in which is a constant calibrated from experimental data and d max 12 is the maximum shear damage. In the following work, a new Eq. (13) is proposed based on the experimental work.

Material
In this study, the Tepex® dynalite 102-RG600(4)/47% material was investigated. This is a woven fabric thermoplastic composite material also known as organo sheet. The fabric is made of roving glass and the matrix material is a PA6 polymer. The weaving of the fabric is in twill form with an equal weight percentage of fibers along the warp and weft, making the longitudinal and transverse mechanical and thermal properties almost equal. (See Fig.1) The material has high specific strength and stiffness along the predefined direction, combined with high toughness and fatigue resistance. It is delivered in the form of sheets and can be thermoformed, making it feasible for mass production and hence rendering it as cost-effective. Combining it with short or long fiber-reinforced plastic of the same type in the form of ribs/force transmission elements offers excellent lightweight potential. Some basic mechanical and thermal material properties are given in Table1.

Material Testing
Measuring shear properties of unidirectional composites as well as those of woven fabricreinforced composite materials are relatively complex in comparison to isotropic materials. There are various standardized test methods to determine the intra-laminar shear properties of composite materials. These methods differ in size and geometry of gauge length, complexity of fixture required, and cost. Many of these methods differ in their ability to produce the pure and uniform shear stress state in the material [5][6][7][8][9][10].
The most suitable intra-laminar shear testing method involves the torsional testing of cylindrical samples [6]. A pure shear stress state can be achieved by twisting the tube sample with fibers oriented at 0° or 90°. This test can produce results with high accuracy. However, the manufacturing of a cylindrical specimen out of organo sheet is very difficult. Therefore, cylindrical specimens are mostly used for testing pipe components in which specimen geometry resembles the actual component.
The v-notch shear test method [5], also known as the Iosipescu Shear test, is performed on a flat, rectangular specimen with centrally located v-notches. It is an asymmetric, fourpoint bending test in which the location of the failure is pre-determined by v-notches. A constant shear force and zero bending moment in the vertical symmetric plane of the specimen are achieved. Thus, a pure shear stress state is achieved in an infinitesimal area between the two notches. However, the shear stress state close to the notches and the middle of the specimen is not uniform. The fundamental assumption of the v-notch shear test method is that material is relatively homogeneous with respect to the size of the specimen. The roving glass fabric has coarser fabric features with a repeating unit that is larger than the distance between the v-notches. Because of this, the Iosipescu test cannot be used for measuring the shear properties of Tepex dynalite 102-RG600 (47%).
The Shear Frame test [7] uses a rectangular specimen with a gauge area of 105 mm × 105 mm. The specimen is clamped in a frame with two identical halves with the help of bolts. Both halves can deform the specimen into a rhombic shape. When loaded in tension, this fixture introduces shear forces in the specimen. Due to a larger gauge area of frame shear, the specimen's out-of-plane bending is observed in specimens with low bending stiffness. The maximum commercially available thickness of the material under consideration in this paper was 4 mm, which does not provide enough bending stability.
Therefore, in this study, three types of shear testing methods were used to test the Tepex dynalite 102-RG600/47%. A schematic illustration of loading type with initial (solid lines) and deformed shapes (dotted lines) is summarized in the second row of Table2. In the Rail shear method [8], the normal stress levels along the fiber direction 1 , 2 are theoretically equal to zero, which results in the pure shear stress-strain response. For the Tensile Shear method [9], in addition to the shear stress tensile or   positive normal, stresses along the fiber direction are also present. In the Compression Shear test method [10,27] shear deformation of the material is accompanied by compression or negative normal stresses along the fiber direction. With the Tensile Shear and Compression Shear methods, the absolute value of normal stresses is equal to shear stress.

Rail Shear Test
The two-rail shear test was carried out according to ASTM D 4255 [8]. This test can be performed on a flat specimen with simple geometrical features, as shown in Table2. The test specimen for this study had six holes of 13 mm diameter each for bolts to pass through them. The specimens were cut with a water jet cutting machine according to the dimensions given in Table2. The test fixture required is shown in Fig.2. The specimen was clamped between the two pairs of rails with bolt and tensile loading applied to the rails, which introduced shear strain in the laminate. Compression loading could have also been applied but larger rails would have been needed to avoid contact between the specimen's upper/lower edges with the fixture at large deformations. The bolts were tightened with a torque of 100 Nm to give a firm grip. The rails of the fixture were 30 mm wide, which covered 60 mm width of the specimen. The visible area of the specimen for optical strain measurement was 14 mm × 150 mm.
The test was done with a constant cross head speed of 3 mm/min. The strain was measured with an optical strain measuring system, Aramis 3D Camera System manufactured by GOM Germany It uses digital image correlation (DIC) technique to measure the strain on the surface of material. Shear stress was calculated as: Here, F is force measured during the test and l is the length of the specimen, which was equal to 150 mm. h is the thickness of the specimen.
In this test method, the upper and lower free edges of the specimen had, theoretically, zero shear stress because no those edges were not constrained. The optically measured strain in Fig. 3 was not uniform. For results evaluation, the average value of the measured strain was calculated. This method assumes a pure shear stress state in the material, although an off-axis load of two rails introduces a small tensile load in the material. In Sect. 4.4, the results are presented and compared with other methods in which multi-stress state loading was applied to the material.

Tensile Shear Test
A tensile shear test was conducted according to ISO 14129 [9]. This test resembles standard tensile testing except that fibers are oriented at ± 45° from the direction of loading. The ease of the test and the requirement of no special fixture makes it the preferred method for measuring the shear properties of composite materials. Using this method, normal tensile stresses are also present in the material in addition to the shear stress during the testing. A complex stress state exists close to the free edges of the specimen. The resulting failure stress is not an exact representation of pure shear response. Although this test method is not recommended for shear response evaluation beyond 5% shear strain due to fiber rotations, with consideration for crash testing it is necessary to continue loading until failure occurs.
The specimen had a gauge length of 150 mm. End tabs measuring 50 mm × 25 mm × 1.5 mm with a fiber orientation of ± 45° angle were glued with fast-curing adhesive Sicomet 77 after roughening the surface. The specimen and test setup are shown in Fig.4.
A stochastic pattern with black and white paint was used for DIC strain measurement. Tests were carried out at the speed of 2 mm/min and strain was measured optically in both the X and Y directions. At the end, shear stress and shear strain were calculated as follows: The diagonal patterns in the measured strain shown in Fig.5 were due to the fibers' direction in ± 45° at the surface. Since strains were only measured on the visible surface, for the test evaluation average value of the strain over the whole surface was calculated.

Compression Shear
There is no standardized test available for compression shear testing. In ISO 14129 [9], the tensile test is performed on a specimen with diagonally oriented fibers to evaluate shear properties. From a mechanic's point of view, the same can also be obtained by performing a compression test on a specimen with diagonally oriented fibers. The  only difference would the change of normal stresses from tensile to compression. The specimen used for this test was 4 mm thick to achieve the maximum possible buckling stability. The dimensions of the test specimen were 146 mm × 25 mm × 4 mm with a gauge length of 18 mm. The combined loading compression (CLC) fixture was used for compression testing of ± 45° laminate. This fixture had four blocks and the specimen was held between the pair of lower and upper blocks with the help of screws (see Fig.6) Four guide pillars were used for the alignment, which were press-fit in the lower blocks and friction-free movement was made possible by the use of a linear ball bearing in the upper blocks. The test specimen was fixed in the fixture with eight screws tightened at 10 Nm with a torque gauge. The upper and lower ends of the specimens were made even with the blocks so that the compression load was applied through the machine head, combined with shear forces due to friction between the specimen and the fixture blocks.
The tests were carried out at a speed of 1 mm/min.
In Fig. 7 the direction of the fibers in diagonal direction can be clearly seen from the measured optical strain pattern. A homogenized strain value was calculated from the average value of the strains.
The compression tests were, in general, difficult to perform in terms of getting less scattered resulting data due to the sensitivity of buckling to the fiber waviness and geometrical imperfection. In axial crash loading, since most of the materials are under compression load, it is necessary to adopt an average curve for compression shear for

Results and Comparison
The shear properties and stress-strain curve results, measured by use of the three testing methods discussed above, are presented in Table3 and Fig.8.
It can be seen that the In-plane shear stress-strain responses were similar for the rail shear method and the tensile shear test. The effect of the multi-axial stress state in the tensile shear testing did not affect the stiffness. However, it led to an increase in the failure strength and strain value. Therefore, an appropriate failure criterion must also be chosen when dealing with complete failure of the material. However, the compressive shear curve was significantly lower than the tensile shear curve due to the difference in fiber rotation. A schematic of the fiber rotation is shown in Fig.9 with solid lines representing the initial position of the fibers and dashed lines representing the fiber orientation after the deformation-driven rotation of the fibers. In tensile loading, the fibers tended to rotate toward the loading direction, which increased the material's ability to bear higher loads. In contrast, under compression loading of angle-ply laminates, the fibers rotated away from the loading direction and hence reduced the load-bearing capacity.

Cyclic Testing
It can be seen from the shear stress-strain plot in Fig.8 that the shear response of the material entered into a non-linear region at a shear strain > 5%. To quantify the plastic deformation and stiffness reduction due to damage, monotonic tests are not sufficient. For this purpose, cyclic tests were carried out with the same geometry and setup. Since the results produced by the rail shear method and the tensile shear test method produced the same results (with the exception of the failure strain and strength for modeling purpose), a full shear stress strain was required, so cyclic tests were performed with tensile shear and compression shear methods only. In cyclic loading, a force of 1000 N was increased in every cycle of loading before unloading. The results are shown in Fig.10 with the dotted lines representing tensile shear compared to compression shear (shown in solid black).
Two important things can be seen directly from the results shown in Fig.10.. The first is that with each loading/unloading cycle, the amount of inelastic strain pl 12 = 2 pl 12 increased; this was plastic deformation. The second is that the average slope of loading and unloading in every cycle decreased with each higher cycle. Ignoring the hysteresis, a simplified illustration is shown in Fig.11.
This stiffness degradation is termed as damage d 12 , which is a phenomenological parameter to quantify the reduction in slope of shear modulus. If the shear elastic modulus in the first cycle was G o 12 , corresponding to undamaged material reduced to G n 12 in the n th cycle, then the shear damage can be calculated as [19,20,24]: The damage is an indirect measurement of the damage area inside the material caused by microcracks. Consequently, the "effective area" inside the material is somewhat less than the original area of the material. This leads to the introduction of the term The effective stress is different from the true stress, which does not account for the damage inside the material. Since only the undamaged part of the material transmits the stress, the constitutive behavior of the material, including the failure criterion and plasticity law, must be determined in terms of effective stress rather than true stress.   In-plane cyclic shear behavior measured by tensile and compression shear methods Therefore, to further compare the difference between tensile shear and compression shear, data acquired from cyclic testing could be used to compare the damage evolution and hardening curves of both tests, as seen in Fig.12 and Fig.13.
In Fig.12, the curves intersecting with the x-axis corresponding to zero shear damage mark the damage onset shear stress S.
In light of the results shown in Fig.12 and Fig.13 it is evident that on the macroscale, the shear curve followed two different hardening curves as well as different damage evolution laws depending upon the normal stress state in the material. A non-linear function must be used to describe the shear damage evolution law and plasticity. The material model should be able to distinguish the shear behavior in the presence of tensile or compression loading.

The Extended CDM Model
As seen in Sect. 4.4, that accumulation of plastic strain and stiffness degradation due to damage were dominating factors in the non-linear shear behavior of the woven fabric composites. It was found that the presence or absence of tensile normal stress or compression normal stress resulted in profound differences in measured shear properties. Therefore, CDM (see Sect. 2) must be able to distinguish between tensile shear and compression shear. A simple method, then, is proposed to detect and choose the materials' constants for shear based on the loading direction: Further, for calibration of shear damage parameters, experimental data were used, as proposed by Johnson [20], according to Eq. (8) in Fig. 14.
It is obvious that the logarithmic function proposed by Johnson [20] does not fit the experimental data well. Instead an exponential function for calibrating the damage is proposed as: in which m, n, and o are material constants measured from the cyclic shear test. This function restricts the maximum value of damage to d max 12 , which was measured experimentally. Thus, for the organo sheet, the damage evolution function requires 3 parameters to be determined from a plot between the damage threshold and the damage.

Results
The user material model code was written as Fortran code and implemented as VUMAT in Abaqus. A pseudo-code for implementation of the algorithm is given in the Appendix. The verification of the results was done with cyclic loading of a single element with the fiber  Fig.15. For better graphical presentation, only two cycles of the simulation results are shown. As can be seen from the results shown in Fig.15, there was reasonable agreement between the simulation results and the experimental data. The plasticity and stiffness degradation were captured very well in the material model.
For further verification of the user material model, monotonic simulation results of the tensile shear and compression shear of the single element model were also compared with the experimental data shown in Fig.16. The output of solution-dependent variables enabled the comparison of inelastic strain and damage. The experimental data were from the cyclic tests described in Sect. 4, from which the material parameters were calibrated.
It was noted that the simulation results agreed with experimental results. Importantly, the use of the user material model enabled detection of the tensile/compression loading by itself; the input material parameters were then used accordingly. As such, the user-defined model can be used for the simulation of further models.

Validation
After verification of the model on single element simulations, simulation results were validated on three-point bending of ± 45° laminate of four plies. The test specimens had rectangular dimensions of 100 mm × 25 mm × 2 mm. The impactor and two support rollers had a 10 mm diameter with supports at a distance-defining span length of 60 mm. A displacement-controlled loading was applied at the speed of 5 mm/min until the force level reached the maximum and then dropped. It was kept even at extreme bending, so that the specimen could not come in contact with the internal part of fixture other than the supports and impactor.
This test was chosen because in three-point bending, the lower half of the material was under tensile loading and the upper half experienced compression stress at the same time. Due to the diagonal orientation of the fibers, the material also underwent shear deformation. Therefore, the lower half of the sample, under bending, experienced tensile shear stress while at the same time the upper half of the sample experienced compression shear stress. Therefore, this simple and comprehensive test was used to verify whether the user material subroutine can simultaneously detect the loading type in finite elements. It also pointed to ways in which the user material model was able to automatically detect the compression and tensile shear loading and how it improved the simulation results.
A simulation model is shown in Fig.17 with a rigid impactor and supports. The impactor was modeled as a cylindrical form as only the radial part of the impactor came in contact with the specimen. To reduce the computational time of explicit simulation, the simulation time was reduced to 20 ms. To avoid the dynamic effects, the displacement boundary conditions were applied on the rigid impactor using Smooth amplitude. This method gradually ramped up the deformation and resulted in the avoidance of undesirable inertial effects. To avoid the slippage of the specimen in the z-direction, the tangential frictional coefficient was set equal to 0.03. To confirm the quasi-static response of the simulation and to avoid unwanted errors, the energy output was requested in the simulation model, as shown in Fig.18. The total strain energy E I was almost equal to external work E w . Energy lost to kinetic energy E KE , frictional dissipation E FD , artificial strain energy E AS , and viscous dissipation E VD were negligible in comparison to external work. Finally, the simulation results were plotted against experimental results shown in Fig.19.
In Fig.19 the simulation curve labelled as direction dependent shows the result where the user material model automatically chose tensile shear parameters for finite elements under tensile loading and compression shear parameters for elements under compression loading. The simulation curve labeled as Tensile Shear shows results where only shear parameters measured from the Tensile Shear Test were used, which was generally used up to now. For deformation up to 7 mm, those test data and simulation data are in good agreement. At higher deformation, the gap between the test data and the simulation data increases but still direction dependent simulations results are closer to the experimental data and better than state-of-the-art. Therefore, it was noted that automatic detection of the loading direction and the selection of corresponding parameters improved the accuracy of the predictions. The need for automatic detection of direction-dependent parameters is more important for the load cases in which off-axis laminates are subjected to compression loading, e.g. with an axial crash of cylindrical tubes.

Conclusions
In this study, the influences of tensile and compression stress in the fiber direction relative to shear behavior were investigated. The rail shear method was used to measure the pure shear behavior of the composite, while tensile and compression testing of ± 45° laminate were used to measure shear behavior in the presence of tensile and compression stresses, respectively. It was found that: • Measured shear properties of woven fabric composite are dependent upon test type; • Shear strength and shear failure strain increase in the presence of tensile normal stress; • Shear resistance of the material significantly reduces in the presence of compression stresses; • The Ramberg-Osgood hardening and exponential function of damage, in terms of effective shear stress, fits well for modeling the shear behavior twill woven fabric composite material; • Mathematically, the sign of trace of strain tensor was used successfully to detect the type of loading, so that appropriate material parameters can be selected for modeling; • The implemented user material model predicted the shear behavior of the material perfectly on the coupon level; and • The simulation of three-point bending of ± 45° laminate proved the effectiveness of the improved CDM model. Fig. 19 Simulation of three-point bending showing improvement in prediction results by using directional-dependent parameters in comparison to parameters measured by tensile shear only small-strain increment at time t = n + 1 . At first, every new increment is supposed to be elastic and the corresponding stress value is termed as Trial Stress ∼ trial n+1 . The trial value is plugged into yield function Eq. (6) to verify that it satisfies the yield condition. If the yield condition is not satisfied, this would mean that the new increment was not a pure elastic step. The return algorithm then calculates the plastic strain − pl n+1 , and the trial stress is corrected (see step 6 in Table4). Table4 summarizes the steps for implementation of the material model as user-defined code.