Compressive Damage and Failure Behavior of Plain-Woven GFRP Sandwich Composite L-Joints of Ships: Experiment and Simulation

Abstract

This paper studies the compressive failure behavior of pw-GFRP (plain-woven glass fiber reinforced polymer) sandwich composite L-joints for ships. Six L-joint specimens were subjected to ultimate compressive tests, which show that damage and stiffness degradation started much earlier than any visible signs. A modified fiber kinking criterion is developed for the compressive failure of pw-GFRP layup structure. Moreover, a hybrid stiffness degradation model is applied to simulate the damage progression. The criterion and the degradation model are implemented in simulation through Abaqus user coding. The simulation results are compared with other existing theories, and the proposed criterion is validated by the similar damage initiation time. The combination of the proposed criterion and hybrid degradation model demonstrates a clear advantage over conventional methods in predicting the compressive damage evolution and ultimate failure behavior of the L-joint structure.

1. Introduction

In recent years, the world’s major naval powers are actively promoting the application of composite structures on warships [1]. These structures are typically larger, more modular, and more platform-oriented than traditional ones. Sandwich composite structures have better bending resistance than conventional laminates. Therefore, they are a common choice for replacing conventional load-bearing structures. L-joint is one of the basic forms of connection in ship structures. Its failure behavior plays a critical role in the response of large composite marine structures under extreme conditions.

Previous studies on the compressive failure behavior of composite L-joints of ship structures are scarce [2]. Qiu Jiabo et al. [3] conducted bending tests on stiffened L-joints. They found that they had high ultimate bearing capacity but also large deformation in the ultimate state. Li et al. [4] analyzed the performance of a new type of L-joint and optimized the design to enhance its bending stiffness by nearly 40%. Zeng et al. [5] studied an L-joint structure similar to the one in this paper. They performed progressive failure simulations of the tensile and compressive ultimate conditions using a modified Hashin criterion combined with an instant stiffness degradation strategy. Their numerical model could predict the failure loads accurately. Shen Wei [6] investigated the same structure as Zeng Haiyan, but focused on the fatigue performance and developed a fatigue life calculation model based on the test results. Most of these studies [2,4,5] concentrated more on ultimate loads, while their simulation results did not match well with the damage process and structural response. It is mainly due to the unique geometry of the L-joint that any influence from the damage at the transitional corner is amplified at the two ends of the joint. Therefore, higher demands are placed on the accuracy and applicability of the damage model.

The main challenge in this study is to develop a suitable macroscopic progressive damage model (PDM) for plain woven glass fiber reinforced polymer (GFRP) under compression loading, which should be able to accurately capture the damage evolution process. The macroscopic PDMs see a single or multiple plies as a whole and consider them as homogeneous anisotropic materials. The primary objective in modelling progressive failure is the adoption of appropriate failure criteria and damage evolution laws. These determine when the material fails and how the structure behaves after failure.

When the research object is a large-scale marine composite structure (such as T-joint, L-joints in this paper, panel-frame structure, etc.), there are many possible failure locations and modes. Therefore, the traditional mode-related criteria are favored by most researchers due to their wide applicability and simplicity. The most famous Hashin failure criteria [7,8,9] emerged and gradually established a dominant position in the grand “World Wide Failure Exercise” (WWFE) [10,11,12], and researchers have also continuously improved and promoted it in specific fields [13,14,15,16]. Shokrieh et al. [13] further extended the Hashin criterion and proposed the Shokrieh-Hashin criterion, which can predict up to seven failure modes in three-dimensional directions. Such methods are simple and easy to implement, yet cannot meet the high accuracy requirements posed by the unique structure at the transitional corner of the L-joints.

Davila and Camanho [17] later developed the LaRC series of failure criteria based on the ideas of Sun, Puck, Hahn, and Johannsson [18]. The LaRC has a large impact due to its relatively simple form and satisfying results. In particular, the LaRC03#6 criterion for fiber compression failure is based on the elastic fiber kinking mechanism. Fiber kinking failure is usually observed in composite materials with high fiber volume fraction under the compressive load [19]. This method evaluates failure by determining whether the matrix shear or tensile fracture occurs near the misaligned fiber bundle, which is a significant improvement over the previous criteria based on maximum stress. Pinho [20] proposed a fiber kinking criterion in three-dimensional space based on this theory, but there are also problems such as the difficulty in determining the initial fiber misalignment angle and the angle of the fiber kinking damage surface; more importantly, the theory is only applicable to unidirectional plates.

In terms of damage evolution laws, existing theories are divided into instantaneous and continuous stiffness degradation models. Instantaneous stiffness degradation decreases the stiffness instantly to a predefined ratio once the corresponding type of damage occurs. It is simple and easy to implement. The downside is that altering element stiffness instantly is a strong nonlinear process, which poses little control over damage evolution behavior. As a result, damage evolution using instantaneous degradation is very dependent on mesh size and local geometry [21,22,23,24]. The continuous stiffness degradation model (CDM) is a method that calculates the energy dissipation based on the damage driving force of the failure criterion. This method was first proposed by Camanho [25] when using cohesive elements to simulate interlayer delamination damage, and then was generalized and applied to various damage evolution processes of composite materials [26,27,28,29,30].

In this paper, the authors first describe and analyze the compressive failure behavior of the L-joint specimens during experimental tests. Then, to better simulate the compressive failure of the L-joints, this paper proposes an applicable compressive fiber kinking criterion for the plain-weave GFRP (pw-GFRP) layup structure. Moreover, a versatile hybrid stiffness degradation model is also proposed to adopt advantages of both instantaneous and continuous models. The combined PDM is implemented in Abaqus through user coding. Finally, the result is compared with both experimental data and simulation results using other conventional methods.

2. Tests

2.1. L-Joint Specifications and Test Setup

The schematic diagram of the L-joint is shown in Figure 1. This joint is a typical form of connection for large-scale composite marine structures. It consists of two reinforced sandwich panels that form a 105° angle at the intersection. The longitudinal stiffener is joined with a rounded transition in the mid corner. The skin of the sandwich structure is made of 3 mm-thick plain-woven GFRP with six layup layers. The core material is PVC (polyvinyl chloride). The detailed material parameters are listed in

Table 1. Material Properties.

Material Properties for pw-GFRP Skin (Material Data are Provided by the Manufacturer, Xianning Haiwei Composite Materials Co., Ltd.)
Elastic Modulus	Fiber Directions $E_1=E_2=18000 \mathrm{M}\mathrm{P}\mathrm{a}$	Layup Direction $E_3=5000 \mathrm{M}\mathrm{P}\mathrm{a}$
Shear Modulus	$\mathrm{Inplain} \mathrm{Direction} G_{12}=3550 \mathrm{M}\mathrm{P}\mathrm{a}$	$\mathrm{Interlaminar} \mathrm{Direction} G_{13}=G_{23}=3550 \mathrm{M}\mathrm{P}\mathrm{a}$
Poisson’s Ratio	$\mathrm{Inplain} \mathrm{Direction} {\nu }_{12}=0.13$	$\mathrm{Interlaminar} \mathrm{Direction} {\nu }_{13}={\nu }_{23}=0.3$
Axile Strength	$\mathrm{Fiber} \mathrm{Compression} X_{\mathrm{f}\mathrm{k}\mathrm{c}}=240 \mathrm{M}\mathrm{P}\mathrm{a}$	$\mathrm{Interlaminar} \mathrm{Compression} Z_{\mathrm{c}}=185 \mathrm{M}\mathrm{P}\mathrm{a}$
Shear Strength	$\mathrm{Inplain} \mathrm{Direction} S_{12}=160 \mathrm{M}\mathrm{P}\mathrm{a}$	$\mathrm{Interlaminar} \mathrm{Direction} S_{13}=S_{23}=80 \mathrm{M}\mathrm{P}\mathrm{a}$
Fracture Energy	$\mathrm{Fiber} \mathrm{Kinking} G^{\mathrm{f}\mathrm{k}\mathrm{c}}=1 \mathrm{J}·{\mathrm{m}\mathrm{m}}^{-2}$	$\mathrm{Interlaminar} \mathrm{Compression} G^{\mathrm{d}\mathrm{c}}=4 \mathrm{J}·{\mathrm{m}\mathrm{m}}^{-2}$
Damage Variable	$\mathrm{Threshold} \mathrm{for} \mathrm{Fiber} \mathrm{Kinking} {DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}=0.6$	$\mathrm{Maximum} \mathrm{Failure} \mathrm{DV} {DV}_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.999$
Material Properties for H100 PVC Foam (Material data are acquired through sample tests.)
Engineering Properties	Elastic Modulus $E=100\mathrm{MPa}$	Poisson’s Ratio $\nu =0.32$
	$\mathrm{Tensile} \mathrm{Yield} \mathrm{Strength} {\sigma }_{\mathrm{t}}=3.10\mathrm{MPa}$
Crushable Foam	Compressive Yield Strength Ratio k	Hydrostatic Compressive Yield Strength Ratio kt
	$k=\frac{{\sigma }_{\mathrm{c}}^0}{P_{\mathrm{c}}^0}\approx 1.625$	$k_{\mathrm{t}}=\frac{P_{\mathrm{t}}}{P_{\mathrm{c}}^0}\approx \frac{{\sigma }_{\mathrm{t}}}{{\sigma }_c^0}=\frac{3.10 \mathrm{M}\mathrm{P}\mathrm{a}}{1.172 \mathrm{M}\mathrm{P}\mathrm{a}}=2.645$
	where ${\sigma }_{\mathrm{c}}^0$ is the initial yield stress in uniaxial compression; $P_{\mathrm{c}}^0$ is the initial yield stress in hydrostatic compression. Hydrostatic compression tests are difficult to conduct. Here k takes the value of 1.625 by referring to this article [33].	where Pt is the yield strength in hydrostatic tension; ${\sigma }_{\mathrm{c}}^0$ = 1.172 MPa is acquired from material tests. The value of Pt and $P_{\mathrm{c}}^0$ are difficult to acquire. Therefore kt is estimated using ${\sigma }_{\mathrm{t}}/{\sigma }_{\mathrm{c}}^0$
Volumetric Hardening	$\mathrm{Initial} \mathrm{yield} \mathrm{stress} {\sigma }_1=1.172 \mathrm{M}\mathrm{P}\mathrm{a}$	$\mathrm{Yield} \mathrm{stress} 2 {\sigma }_2=1.800 \mathrm{M}\mathrm{P}\mathrm{a}$
	$\mathrm{Initial} \mathrm{uniaxial} \mathrm{plastic} \mathrm{strain} {\epsilon }_1=0$	$\mathrm{Uniaxial} \mathrm{plastic} \mathrm{strain} 2 {\epsilon }_2=0.00335$

. The tensile damage behavior and the simulation methods of this joint have been thoroughly discussed in two previous articles [31,32]. This paper mainly focuses on the compressive damage behavior and its simulation methods.

There are two types of specimens with varying corner radii of the transitional area (R). Three specimens have R = 45 mm, numbered CR45-1 to CR45-3, and three specimens have R = 90 mm, numbered CR90-1 to CR90-3. During the L-joint test, the fixture setup clamps the two ends while ensuring its free rotation around the corresponding axis shown in Figure 2a. A mandatory compressive displacement at a speed of 2 mm/min is applied to the upper end until the structure breaks.

2.2. Test Results

Compressive L-joints typically experience failure at the middle strengthener, with initial damage observed on the edges of the strengthener. However, the specific location of the initial damage varies depending on the radius of the transitional area (R).

In the case of CR45 specimens, whitened spots are first observed on the two edges near the midpoint of the transitional arc, as shown in Figure 3b. (Figure 3a presents the original state for comparison). As the load increases, these spots grow into two distinct cracks that expand towards each other on the upper surface of the strengthener (Figure 3c). The L-joint fails quickly once the cracks join into one (Figure 3e) and then extend down both sidewalls of the strengthener from their initial spot (Figure 3d,f,g).

The failure process of CR90 specimens is similar. However, damage starts on both edges but not on the horizontal symmetry plane of the transitional arc (see Figure 4b), Figure 4a shows the original state for comparison). When the two whitened spots are on different sides of the symmetry plane, the extended cracks join to form an “S”-shape, as illustrated in Figure 4b,c,f. When they are on the same side, the final crack is nearly a straight line, as depicted in Figure 4e. Eventually, the crack also propagates along the side walls (Figure 4d,g).

2.3. Test Data Analysis

In addition to the displacement and load data obtained from the actuator of the test machine at the upper end of the L-joints, we also measured the surface strain at several locations on the upper surface of the stiffener along its longitudinal direction. Figure 1 shows that point B is on the horizontal symmetry plane; points A and C are symmetrical to each other (we excluded the strain results at point C from later analysis because they were similar to point A); point D is at the midpoint of the curve length between A and B. The strain data were collected and recorded at a frequency of 5 Hz, consistent with the test machine data. Figure 5 plots the L-D (load-displacement) curve and L-E (load-strain) curves of each specimen.

The manufacturing process disperses the structural properties of resin-based glass fiber composites. Therefore, the L-E curves of the same type of specimens will show significant numerical differences if compared directly. The main concern of this test is the damage process of the specimen. The overall pattern reveals that the L-D and L-E curves are linear when the load force is less than half of the ultimate strength. As the load increases, gradual changes in the corresponding trends indicate that some structural damage has occurred.

For specimen CR45-2, for example, strain B increases more rapidly while strain D decreases after the load force surpasses slightly over 30 kN. The strain trend changes gradually and smoothly in contrast to the sudden jumps under tension [32]. This indicates that the specimen is progressively damaged, but it is hard to pinpoint the exact moment or location of damage initiation. However, the video recording of each test shows the damage on the upper surface of the stiffener. The appearance of “white spots” on this surface coincides with the strain trend alteration in Figure 5. Thus, the load forces corresponding to these spots are recorded as the initial signs of damage during the tests.

Table 2. Test Load Results Summarization (Unit: kN).

	Specimen 1	Specimen 2	Specimen 3	Mean	Std. Deviation
CR45 Damage Observed	30.28	33.74	31.15	31.72	1.800
CR45 Ultimate Strength	46.45	45.82	43.2	45.16	1.724
CR90 Damage Observed	51.8	51.07	56.36	53.08	2.867
CR90 Ultimate Strength	60.88	63.62	61.98	62.16	1.379

lists the load of initial damage observation and ultimate strength for each specimen.

3. Damage Criterion

3.1. Fiber Kinking

The observed failure was directly caused by longitudinal compressive damage of the skin. Recent technological advances have revealed that most fiber compressive failure involves fiber microbuckling [34] and fiber kinking band [35,36], which results from damage of the surrounding matrix that supports the fiber bundles. Some researchers [35,37] argue that initial microscopic defects, such as fiber misalignment, matrix cracking, or interface cracking, exist inherently. Camanho and Davila [17] define these initial defects by a parameter called initial fiber misalignment angle, which is a constant of the material property. The fiber misalignment angle is used in the local stress analysis of the fiber kinking area. It will evolve when the material is under a compressive load on fiber direction.

Figure 6a shows a planar case of unidirectional material subjected to axial compression, with $\theta$ representing the fiber misalignment angle. The fiber direction varies in the local misalignment area. Based on the LaRC03#6 criteria developed by Davila [17], fiber kinking failure is assessed by examining the matrix tensile or compressive damage in this area:

$$\{\begin{array}{c}1\mathrm{st} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} (\mathrm{matrix} \mathrm{tension}) \\ {\Phi }_1^{\mathrm{fkt}}={\left(\frac{{\sigma }_2^{\mathrm{fm}}}{Y_{\mathrm{mt}}}\right)}^2+{\left(\frac{{\tau }_{12}^{\mathrm{fm}}}{S_{\mathrm{L}}}\right)}^2\ge 1,\left({\sigma }_2^{\mathrm{fm}}\ge 0\right)\\ 1\mathrm{st} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} (\mathrm{matrix} \mathrm{compression}) \\ {\Phi }_1^{\mathrm{fkc}}={\left(\frac{\langle \left|{\tau }_{12}^{\mathrm{fm}}\right|+{\mu }_{\mathrm{L}}{\sigma }_2^{\mathrm{fm}}\rangle }{S_{\mathrm{L}}}\right)}^2\ge 1,\left({\sigma }_2^{\mathrm{fm}}\le 0\right)\end{array}$$

(1)

where $S_{\mathrm{T}},S_{\mathrm{L}}$ are the transverse and longitudinal shear strength.

${\sigma }_1^{\mathrm{f}\mathrm{m}}$, ${\sigma }_2^{\mathrm{f}\mathrm{m}}$, and ${\tau }_{12}^{\mathrm{f}\mathrm{m}}$ are stress components on the local misalignment area, which can be calculated from:

$$\begin{array}{l}\{\begin{array}{l}{\sigma }_1^{\mathrm{fm}}=\frac{{\sigma }_1+{\sigma }_2}{2}+\frac{{\sigma }_1-{\sigma }_2}{2}\cos (2\theta )+{\tau }_{12}\sin (2\theta )\hfill \\ {\sigma }_2^{\mathrm{fm}}={\sigma }_1+{\sigma }_2-{\sigma }_1^{\mathrm{fm}}\hfill \\ {\tau }_{12}^{\mathrm{fm}}=-\frac{{\sigma }_1-{\sigma }_2}{2}\sin (2\theta )+{\tau }_{12}\cos (2\theta )\hfill \end{array}\hfill \end{array}$$

(2)

$$\{\begin{array}{l}\theta =\frac{{\tau }_{12}+G_{12}\frac{{\tau }_{12}}{\left|{\tau }_{12}\right|}{\theta }_{\mathrm{i}}}{G_{12}+{\sigma }_1-{\sigma }_2}\\ {\theta }_{\mathrm{i}}={\theta }_{\mathrm{c}}-{\theta }_{\mathrm{r},\mathrm{c}}={\theta }_{\mathrm{c}}\left(1-\frac{X_{\mathrm{fc}}}{G_{12}}\right)\\ {\theta }_{\mathrm{c}}={\tan }^{-1}\left(\frac{1-\sqrt{1-4\left(\frac{S_{\mathrm{L}}}{X_{\mathrm{fc}}}+{\mu }_{\mathrm{L}}\right)\frac{S_{\mathrm{L}}}{X_{\mathrm{fc}}}}}{2\left(\frac{S_{\mathrm{L}}}{X_{\mathrm{fc}}}+{\mu }_{\mathrm{L}}\right)}\right)\end{array}$$

(3)

in which $\theta$ is the fiber misalignment angle, ${\theta }_{\mathrm{i}}$ is the initial fiber misalignment angle, and ${\theta }_{\mathrm{c}}$ is the critical fiber misalignment angle when fiber kinking damage happens. ${\mu }_{\mathrm{L}}$ and ${\mu }_{\mathrm{T}}$ are the static friction coefficients on the assumed matrix cracking surface.

$$\{\begin{array}{l}{\mu }_{\mathrm{T}}=\frac{-1}{\tan \left(2{\phi }_{\mathrm{mc}}\right)}\hfill \\ \frac{{\mu }_{\mathrm{T}}}{S_{\mathrm{T}}}=\frac{{\mu }_{\mathrm{L}}}{S_{\mathrm{L}}}\hfill \end{array}$$

(4)

${\phi }_{\mathrm{m}\mathrm{c}}$ is the fracture angle acquired from the pure transverse loading test, which is experimentally observed to be 53 ± 2° for most composite materials [19,38,39,40]. Here ${\phi }_{\mathrm{m}\mathrm{c}}$ is taken as 53°.

The LaRC03#6 criteria discuss fiber kinking only in terms of planar orientation and assume a vertical matrix fracture plane. This is a reasonable and convenient assumption for common laminate cases where the fiber kinking plane is physically restricted to a single layer. However, for other cases such as unidirectional materials under a complex stress state, the orientations of both the fiber kinking plane and the matrix fracture plane remain unknown.

Pinho [20] proposed a 3D kinking model for unidirectional materials using similar methods.

$$\{\begin{array}{c}\mathrm{First} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} (\mathrm{matrix} \mathrm{tension})\\ {\Phi }_1^{\mathrm{fkt}}={\left(\frac{{\sigma }_{\mathrm{n}}^{\mathrm{fm}}}{Y_{\mathrm{mt}}}\right)}^2+{\left(\frac{{\tau }_{\mathrm{T}}^{\mathrm{fm}}}{S_{\mathrm{T}}}\right)}^2+{\left(\frac{{\tau }_{\mathrm{L}}^{\mathrm{fm}}}{S_{\mathrm{L}}}\right)}^2\ge 1,\left({\sigma }_{\mathrm{n}}^{\mathrm{fm}}\ge 0\right)\\ 1\mathrm{st} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} (\mathrm{matrix} \mathrm{compression})\\ {\Phi }_1^{\mathrm{fkc}}={\left(\frac{\langle \left|{\tau }_{\mathrm{T}}^{\mathrm{fm}}\right|+{\mu }_{\mathrm{T}}{\sigma }_{\mathrm{n}}^{\mathrm{fm}}\rangle }{S_{\mathrm{T}}}\right)}^2+{\left(\frac{\langle \left|{\tau }_{\mathrm{L}}^{\mathrm{fm}}\right|+{\mu }_{\mathrm{L}}{\sigma }_{\mathrm{n}}^{\mathrm{fm}}\rangle }{S_{\mathrm{L}}}\right)}^2\ge 1,\left({\sigma }_{\mathrm{n}}^{\mathrm{fm}}\le 0\right)\end{array}$$

(5)

The stress components of the matrix cracking plane can be calculated as follows:

$$\{\begin{array}{l}{\sigma }_{\mathrm{n}}^{\mathrm{fm}}=\frac{{\sigma }_2^{\mathrm{fm}}+{\sigma }_3^{\mathsf{\alpha }}}{2}+\frac{{\sigma }_2^{\mathrm{fm}}-{\sigma }_3^{\mathsf{\alpha }}}{2}\cos \left(2{\phi }^{\mathrm{fm}}\right)+{\tau }_{23}^{\mathrm{fm}}\sin \left(2{\phi }^{\mathrm{fm}}\right)\\ {\tau }_{\mathrm{T}}^{\mathrm{fm}}=-\frac{{\sigma }_2^{\mathrm{fm}}-{\sigma }_3^{\mathsf{\alpha }}}{2}\sin \left(2{\phi }^{\mathrm{fm}}\right)+{\tau }_{23}^{\mathrm{fm}}\cos \left(2{\phi }^{\mathrm{fm}}\right)\\ {\tau }_{\mathrm{L}}^{\mathrm{fm}}={\tau }_{12}^{\mathrm{fm}}\cos {\phi }^{\mathrm{fm}}+{\tau }_{13}^{\mathrm{fm}}\sin {\phi }^{\mathrm{fm}}\end{array}$$

(6)

The stress components of the fiber misalignment frame can be calculated as follows:

$$\{\begin{array}{l}{\sigma }_1^{\mathrm{fm}}=\frac{{\sigma }_1+{\sigma }_2^{\mathsf{\alpha }}}{2}+\frac{{\sigma }_1-{\sigma }_2^{\mathsf{\alpha }}}{2}\cos \left(2\theta \right)+{\tau }_{12}\sin \left(2\theta \right)\\ {\sigma }_2^{\mathrm{fm}}={\sigma }_1+{\sigma }_2^{\mathsf{\alpha }}-{\sigma }_1^{\mathrm{fm}}\\ {\tau }_{12}^{\mathrm{fm}}=-\frac{{\sigma }_1-{\sigma }_2^{\mathsf{\alpha }}}{2}\sin \left(2\theta \right)+{\tau }_{12}^{\mathsf{\alpha }}\cos \left(2\theta \right)\\ {\tau }_{23}^{\mathrm{fm}}={\tau }_{23}^{\mathsf{\alpha }}\cos \theta -{\tau }_{13}^{\mathsf{\alpha }}\sin \theta \\ {\tau }_{13}^{\mathrm{fm}}={\tau }_{13}^{\mathsf{\alpha }}\cos \theta \end{array}$$

(7)

Since the initial fiber misalignment angle reflects the initial defect of the material, in this 3D scenario, ${\theta }_{\mathrm{i}}$ may exist in any direction within the cone shown in Figure 6b. The most favorable situation for fiber kinking is when ${\theta }_{\mathrm{i}}$ is aligned with the direction of fiber rotation under the current stress state. Therefore, fiber kinking is most likely to occur on the plane where ${\sigma }_2^{\mathsf{\alpha }}$ reaches the maximum principal stress in (2, 3) plane. Then:

$$\alpha =\frac{1}{2}{\tan }^{-1}\frac{2{\tau }_{23}}{{\sigma }_2-{\sigma }_3}$$

(8)

3.2. Fiber Kinking Damage Criterion for Plain-Woven Fabrics

Many current macroscopic PDM methods for pw-FRP materials adopt the Max-Stress criterion [5] to deal with compressive damage along the fiber direction. However, some researchers have observed that fiber kinking can occur across multiple layers and cause macroscopic compression failure in pw-FRP, even though the fibers are orthogonally woven. Therefore, the Max-Stress criterion may not be adequate for predicting this failure mode. Pinho et al. [20] and Bishara et al. [41] were the first to propose this theory based on theoretical analysis. Opelt et al. [42] confirmed this hypothesis by directly observing the kinking band in a plain-woven carbon fiber reinforced polymer (pw-CFRP). Toribio et al. [43] also reported that fiber kinking occurred during the compressive test of notched pw-CFRP materials. In his research, a photo taken under microscope clearly shows the kinking band of the warp fiber bundle cuts through the weft one. Khan et al. [44] further investigated this phenomenon and summarized its damage and failure mechanism. Therefore, it is reasonable to construct a failure criterion for fiber compression following the principle of fiber kinking.

In the 3D kinking model for unidirectional materials, the angle ${\theta }_{\mathrm{i}}$ depends on the local stress state to account for the most favorable scenario for fiber rotation. However, this is not the case for pw-FRP. Figure 7 shows an illustration of the pw-FRP structure, where the fiber bundles are intertwined due to its unique stacking technique. This intertwined geometry causes a significant initial fiber misalignment along the layup direction (on the (1, 3) or (2, 3) plane depending on the fiber direction, see Figure 7). This effectively restricts the fiber misalignment angle on the same direction, which means, for pw-FRP, the angle $\alpha =90°$ in Figure 6b. Therefore, the following are the 3D fiber kinking criteria for pw-FRP:

$$\{\begin{array}{c}\mathrm{First} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} \mathrm{for} \mathrm{pw}\text{-}\mathrm{FRP} (\mathrm{matrix} \mathrm{tension})\\ {\Phi }_1^{\mathrm{fkt}}={\left(\frac{{\sigma }_2^{\mathrm{fm}}}{Z_{\mathrm{t}}}\right)}^2+{\left(\frac{{\tau }_{12}^{\mathrm{fm}}}{S_{13}}\right)}^2+{\left(\frac{{\tau }_{23}^{\mathrm{fm}}}{S_{23}}\right)}^2\ge 1,\left({\sigma }_2^{\mathrm{fm}}\ge 0\right)\\ \mathrm{First} \mathrm{direction} \mathrm{fiber} \mathrm{kinking} \mathrm{for} \mathrm{pw}\text{-}\mathrm{FRP} (\mathrm{matrix} \mathrm{compression})\\ {\Phi }_1^{\mathrm{fkc}}={\left(\frac{\langle \left|{\tau }_{12}^{\mathrm{fm}}\right|+{\mu }_{13}{\sigma }_2^{\mathrm{fm}}\rangle }{S_{13}}\right)}^2+{\left(\frac{\langle \left|{\tau }_{23}^{\mathrm{fm}}\right|+{\mu }_{23}{\sigma }_2^{\mathrm{fm}}\rangle }{S_{23}}\right)}^2\ge 1,\left({\sigma }_2^{\mathrm{fm}}\le 0\right)\end{array}$$

(9)

where:

$$\{\begin{array}{l}{\sigma }_1^{\mathrm{fm}}=\frac{{\sigma }_1+{\sigma }_3}{2}+\frac{{\sigma }_1-{\sigma }_3}{2}\cos \left(2\theta \right)+{\tau }_{13}\sin \left(2\theta \right)\\ {\sigma }_2^{\mathrm{fm}}={\sigma }_1+{\sigma }_3-{\sigma }_1^{\mathrm{fm}}\\ {\tau }_{12}^{\mathrm{fm}}=-\frac{{\sigma }_1-{\sigma }_3}{2}\sin \left(2\theta \right)+{\tau }_{13}\cos \left(2\theta \right)\\ {\tau }_{23}^{\mathrm{fm}}={\tau }_{23}\cos \theta -{\tau }_{12}\sin \theta \\ {\tau }_{13}^{\mathrm{fm}}=-{\tau }_{12}\cos \theta \end{array}$$

(10)

in which, fiber misalignment angle can be calculated using Equations (3) and (4).

4. Material Degradation Model

4.1. Damage Variables and the Damaged Stiffness Matrix

Continuum damage mechanics (CDM) [45] is adopted to simulate the damage progression. When a certain failure criterion is met and damage appears, the extent of material stiffness reduction in the relevant direction is described by the damage variable (${DV}_{\mathrm{i}}$, the subscript indicates material direction). It is generally believed that when damage occurs in a certain direction, it will also cause the shear stiffness in the relevant direction to decrease. For instance, the shear stiffness in (1, 2) plane should be affected by both ${DV}_1$ and ${DV}_2$. However, using $(1-{DV}_1)·(1-{DV}_2)$ directly to describe the shear stiffness loss is obviously too aggressive. In this case, the square root of this value ($\sqrt{(1-{DV}_1)(1-{DV}_2)}$) seems more reasonable. Then, the damage stiffness matrix ${\mathit{C}}^{\mathrm{D}\mathrm{m}\mathrm{g}}$ is constructed as:

$${\mathit{C}}^{\mathrm{D}\mathrm{m}\mathrm{g}}=\left[\begin{array}{cccccc}{k}_1{C}_{11}& \sqrt{k_1{k}_2}{C}_{12}& \sqrt{k_1{k}_3}{C}_{13}& 0& 0& 0\\ & k_2{C}_{22}& \sqrt{k_2{k}_3}{C}_{23}& 0& 0& 0\\ & & k_3{C}_{33}& 0& 0& 0\\ & & & k_{12}{C}_{44}& 0& 0\\ & \mathrm{s}\mathrm{y}\mathrm{m}& & & k_{23}{C}_{55}& 0\\ & & & & & k_{13}{C}_{66}\end{array}\right]$$

(11)

where $C_{ij}$ are the elements of the undamaged stiffness matrix; $k_1=(1-{DV}_1)$, $k_2=(1-{DV}_2)$, $k_3=(1-{DV}_3)$; $k_{12}=2k_1{k}_2/(k_1+k_2)$, $k_{23}=2k_2{k}_3/(k_2+k_3)$, $k_{13}=2k_1{k}_3/(k_1+k_3)$.

The variable $k_{\mathrm{i}}$ here is often referred to as the residual stiffness ratio. Using $k_{\mathrm{i}}$ instead of ${DV}_{\mathrm{i}}$ in ${\mathit{C}}^{\mathrm{D}\mathrm{m}\mathrm{g}}$ can result in a simpler form, but the concept of ${DV}_{\mathrm{i}}$ is more commonly used by researchers.

4.2. Damage Evolution

With the form of ${\mathit{C}}^{\mathrm{D}\mathrm{m}\mathrm{g}}$ established, the degree of stiffness loss can be measured by ${DV}_{\mathrm{i}}$. When damage occurs, the value of ${DV}_{\mathrm{i}}$ follows a stiffness degradation law. In continuum degradation, $DV$ varies from 0 to 1 continuously as the local fracture energy is absorbed. The most common continuum strategy is linear degradation shown in Figure 8. This model was proposed by Camanho et al. [46] for the simulation of interlayer delamination, which can produce good results for materials with significant plastic characteristics. However, for many fiber-reinforced materials, their fracture behavior before failure is not entirely “continuum”. The fiber kinking band, which is often observed after compressive failure, is generally considered to be formed when the structure collapses suddenly under compression.

Therefore, a hybrid degradation model is adopted to address this problem. The main idea is to use a continuum model after damage initiation until $DV$ reaches a critical threshold ${DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}$. Then, it is assumed that the material cannot sustain any more load and instant degradation occurs to change the $DV$ to 1 or another preset constant depending on the final failure state.

When linear weakening is adopted in the continuum part of the hybrid model, the relationship between equivalent stress and displacement is shown in Figure 9.

In this model, when $X_{\mathrm{e}\mathrm{q}}=X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}={\epsilon }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}·L$, it is the critical point of damage initiation; then,

$$\{\begin{array}{cc}\hfill D{V}^{\mathrm{i}}& =0\hfill \\ \hfill E& ={\sigma }_{\mathrm{eq}}^{\mathrm{i}}/{\epsilon }_{\mathrm{eq}}^{\mathrm{i}}={\sigma }_{\mathrm{eq}}^{\mathrm{i}}\cdot L/X_{\mathrm{eq}}^{\mathrm{i}}\hfill \end{array}$$

(12)

When $X_{\mathrm{e}\mathrm{q}}=X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}={\epsilon }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}·L$, it is the critical point before complete material failure; then

$$\{\begin{array}{cc}\hfill D{V}^{\mathrm{f}}& =D{V}_{\mathrm{val}}\hfill \\ \hfill E^{\mathrm{f}}& ={\sigma }_{\mathrm{eq}}^{\mathrm{f}}/{\epsilon }_{\mathrm{eq}}^{\mathrm{f}}={\sigma }_{\mathrm{eq}}^{\mathrm{f}}\cdot L/X_{\mathrm{eq}}^{\mathrm{f}}\hfill \end{array}$$

(13)

based on the definition of $DV$, $E^{\mathrm{f}}=E·\left(1-{DV}^{\mathrm{f}}\right)$. Therefore,

$$\frac{{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}}{X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}}=\frac{{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}}{X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}}\left(1-{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}\right)$$

(14)

Since fracture energy $G_c$ equals the enclosed shaded area under the ${\sigma }_{\mathrm{e}\mathrm{q}}-X_{\mathrm{e}\mathrm{q}}$ curve in Figure 9, then

$$G_c=\frac{1}{2}{X}_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}+\frac{1}{2}\left({\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}+{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}\right)\left(X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}\right)$$

(15)

Combine Equations (14) and (15):

$$\{\begin{array}{l}{X}_{\mathrm{eq}}^{\mathrm{i}}\cdot {\sigma }_{\mathrm{eq}}^{\mathrm{f}}=\frac{\sqrt{{\left(D{V}_{\mathrm{val}}\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}\right)}^2+8\left(1-D{V}_{\mathrm{val}}\right)\cdot G_c\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}}-D{V}_{\mathrm{val}}\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}}{2}\\ X_{\mathrm{eq}}^{\mathrm{f}}\cdot {\sigma }_{\mathrm{eq}}^{\mathrm{i}}=\frac{\sqrt{{\left(D{V}_{\mathrm{val}}\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}\right)}^2+8\left(1-D{V}_{\mathrm{val}}\right)\cdot G_c\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}}-D{V}_{\mathrm{val}}\cdot X_{\mathrm{eq}}^{\mathrm{i}}{\sigma }_{\mathrm{eq}}^{\mathrm{i}}}{2\left(1-D{V}_{\mathrm{val}}\right)} \end{array}$$

(16)

For any point between $X_{\mathrm{e}\mathrm{q}}=X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}$ and $X_{\mathrm{e}\mathrm{q}}=X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}$, $E^{\prime }={\sigma }_{\mathrm{e}\mathrm{q}}/{\epsilon }_{\mathrm{e}\mathrm{q}}$; since $DV=1-E^{\prime }/E$, then

$$DV=\frac{X_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}(X_{\mathrm{e}\mathrm{q}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}})}{X_{\mathrm{e}\mathrm{q}}(X_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}})}$$

(17)

where $X_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}$ is the interception of the linear degradation line on the $X_{\mathrm{e}\mathrm{q}}$ axis. So,

$$X_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}=\frac{X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}}{{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}-{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}}$$

(18)

Substitute $X_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}$ into Equation (17), then

$$DV=\frac{\left(X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}\right)}{\left(X_{\mathrm{e}\mathrm{q}}^{\mathrm{f}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}-X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}{\sigma }_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}\right)}·\left(1-\frac{X_{\mathrm{e}\mathrm{q}}^{\mathrm{i}}}{X_{\mathrm{e}\mathrm{q}}}\right)$$

(19)

Using the conversion method brought up by Kang et al. [29],

$$\{\begin{array}{l}{X}_{\mathrm{eq}}^{\mathrm{i}}=X_{\mathrm{eq}}/\sqrt{\Phi }\\ {\sigma }_{\mathrm{eq}}^{\mathrm{i}}={\sigma }_{\mathrm{eq}}/\sqrt{\Phi }\end{array}$$

(20)

then

$$DV=\frac{{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}·\left(1-1/\sqrt{\Phi }\right)}{1-\frac{2·\left(1-{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}\right)}{\sqrt{{{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}}^2+8·\left(1-{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}\right)·G_c\Phi /\left(X_{\mathrm{e}\mathrm{q}}{\sigma }_{\mathrm{e}\mathrm{q}}\right)}-{DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}}}$$

(21)

This expression calculates the damage variable of hybrid stiffness degradation when linear weakening is used in the continuum stage. It should be noted that if the calculated result is smaller than threshold ${DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}$, then $DV=0$.

The form of $X_{\mathrm{e}\mathrm{q}}{\sigma }_{\mathrm{e}\mathrm{q}}$ depends on the form of the failure criterion $\Phi$. For a quadratic stress form such as

$${\Phi }_{\mathrm{i}}^{\mathrm{M}}={\left({\sigma }_{\mathrm{i}}/A\right)}^2+{\left({\tau }_{\mathrm{i}\mathrm{j}}/S_{\mathrm{i}\mathrm{j}}\right)}^2+{\left({\tau }_{\mathrm{i}\mathrm{k}}/S_{\mathrm{i}\mathrm{k}}\right)}^2$$

(22)

Fang [30] provided the following calculation method:

$$X_{\mathrm{e}\mathrm{q},\mathrm{i}}·{\sigma }_{\mathrm{e}\mathrm{q},\mathrm{i}}={\sigma }_{\mathrm{i}}{\epsilon }_{\mathrm{i}}+{\tau }_{\mathrm{i}\mathrm{j}}{\epsilon }_{\mathrm{i}\mathrm{j}}+{\tau }_{\mathrm{i}\mathrm{k}}{\epsilon }_{\mathrm{i}\mathrm{k}}$$

(23)

For shear stress criterion such as

$${\Phi }_{\mathrm{i}}^{\mathrm{M}}={\left(⟨\left|{\tau }_{\mathrm{i}\mathrm{j}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{L}}{\sigma }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩/S_{\mathrm{i}\mathrm{j}}\right)}^2+{\left(⟨\left|{\tau }_{\mathrm{j}\mathrm{k}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{T}}{\sigma }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩/S_{\mathrm{j}\mathrm{k}}\right)}^2$$

(24)

it can also be constructed following the same form as

$$X_{\mathrm{e}\mathrm{q},\mathrm{i}}·{\sigma }_{\mathrm{e}\mathrm{q},\mathrm{i}}=⟨\left|{\epsilon }_{\mathrm{i}\mathrm{j}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{L}}{\epsilon }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩·⟨\left|{\tau }_{\mathrm{i}\mathrm{j}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{L}}{\sigma }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩+⟨\left|{\tau }_{\mathrm{i}\mathrm{j}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{T}}{\sigma }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩·⟨\left|{\epsilon }_{\mathrm{j}\mathrm{k}}^{\mathrm{f}\mathrm{m}}\right|+{\mu }_{\mathrm{T}}{\epsilon }_{\mathrm{j}}^{\mathrm{f}\mathrm{m}}⟩$$

(25)

5. Simulation Setup

The quasistatic approach is used to perform the simulations with the “Explicit Module” in FE software Abaqus. The explicit algorithm is advantageous for dealing with highly nonlinear calculations. It avoids convergence issues and always produces a result, while giving an accurate estimate of time consumption after a few initial steps.

The material model based on the previously described theories is incorporated into the FE code using user defined subroutine VUMAT. The material parameters for the skin and PVC core are shown in Table 1. The additional parameter ${DV}_{\mathrm{v}\mathrm{a}\mathrm{l}}$ is set to 0.6 based on the authors’ experience. Max-stress criteria with linear stiffness degradation is used to simulate potential delamination failure. The bonding layer is neglected in the simulation due to the compressive stress state and the absence of noticeable debonding damage during the test. Crushable foam with volumetric hardening is applied for the H100 PVC Core.

The FE model (see Figure 2b,c) did not include the transverse strengthener at the L-joint ends, and the end sections’ displacement was kinetically coupled with the center point of the respective rotational axis. The actual specimens’ edges were not perfectly square with sharp corners. This is because the woven cloth skin created a local arc transition at the sharp edges, resembling rounded corners, as shown in Figure 10a. This feature was important because simulations with sharp corners showed abnormal instability failure along the edges as in Figure 10b, which was never seen in the test. Therefore, the FE model used rounded corners with $r=5\mathrm{m}\mathrm{m}$ (average of 10 measurements). Additionally, sensitive analysis on the mesh size has been conducted, and the appropriated mesh length was settled on 1 mm on the upper stiffener to achieve adequate accuracy and decent calculation speed.

6. Results and Analysis

6.1. Comparison between Simulation and Test Results

The simulated failure process of CR45 and CR90 joints is presented in Figure 11 and Figure 12. The first observed damage spots correspond with the simulated high $DV$ areas. The skin cracks that lead to final failure generally match.

The simulated failure process of CR45 joints is shown in Figure 11. Before the final failure, fiber kinking damage occurs (Figure 11(a.1)) at the edge of the upper surface of the stiffener and spreads (Figure 11(a.2 and a.3)). A large area of skin whitening (Figure 11(b.2)) is also observed at this location when comparing the test photos of the damage appearance and expansion with the undamaged state (Figure 11(b.1)). In the simulation results, the area where complete fiber kinking failure occurs (plotted in red in Figure 11(a.4), $DV=1.0$) has a central crack on the upper surface of the stiffener, which extends to its sidewall. Similar cracks (Figure 11(b.3–b.5)) were found upon structural failure during the test.

Figure 12 shows the simulated failure process of CR90 joints. Compared with the CR45 joints, the fiber kinking damage areas with full propagation (Figure 12(a.3)) are much longer, and the initial damaged areas with $DV>0.3$ in Figure 12(a.3) match the whitened areas in the damage extension test photo (Figure 12(b.2)). The simulation results in two essentially symmetrical cracks. The initial locations of these cracks agree with those observed in the test, but they do not form an S-shaped intersection (Figure 12(b.3)) or a single-sided straight-through crack (Figure 12(b.4)). Rather, they create two symmetrical straight-through cracks (Figure 12(a.4)). This discrepancy may be attributed to the perfect symmetry of the finite element model.

Figure 13 provides a direct comparison between the test and simulated L-D curve results (between “tests” and “Our Method” in Figure 13). The marked load force data show good agreement in terms of final failure. The characteristic of the changing trend is consistent with the tests, which is continuous stiffness loss until final failure. Figure 5 shows the simulation results of displacement, load force, and strain data, in which the corresponding force of fiber kinking damage initiation is marked. The corresponding force for damage observed in tests is always higher than that of fiber kinking damage initiation in simulations. This is expected since damage can only be observed after it has propagated to a certain point.

6.2. Simulation Results Comparison between Theories

The most common method of simulating compressive failure along the fiber direction for pw-FRP structures is to apply the maximum stress criterion [2,3,4,6]. However, some researchers have used the quadratic stress criterion to account for the effect of transverse shear stress on fiber compression [5,13].

Table 3. Compressive criteria included in the comparison.

Damage Criteria	Governing Equation
Max Stress/Hashin Criteria [8] Fiber Compression	${\left(\frac{{\sigma }_1}{X_{\mathrm{c}\mathrm{m}}}\right)}^2\ge 1,({\sigma }_1<0)$
Shokreih-Hashin Criteria [13] Fiber Compression	${\left(\frac{{\sigma }_1}{X_{\mathrm{c}\mathrm{m}}}\right)}^2+{\left(\frac{{\tau }_{12}}{S_{12}}\right)}^2+{\left(\frac{{\tau }_{13}}{S_{13}}\right)}^2\ge 1,({\sigma }_1<0)$
Fiber Kinking for pw-FRP	Equation (9)

lists the criteria used in the comparison, and each one is paired with instant, linear, or hybrid stiffness degradation.

The results of the comparison are shown in Figure 13. The damage initiation load for each criterion is quite close (within 4% of the corresponding load of observed damage), and most of the difference lies in the damage evolution process after damage appears. For any criterion paired with instant degradation, the cross-section of the stiffener collapses instantly after the first sign of damage, resulting in an early and abrupt failure, which clearly differs from the tests.

Simulations using Max-Stress or Shokreih-Hashin criteria have unsatisfactory results. Although hybrid degradation shortens the damage evolution process compared with linear degradation, this effect is not obvious for the CR45 model. Moreover, all the resulting L-D curves fluctuate violently for a long time before final failure, which makes it impossible to obtain a stable and reliable ultimate load.

The proposed fiber kinking criterion for pw-FRP eliminates the fluctuation that occurs when using other criteria. The L-D curve with linear degradation shows a longer damage evolution process than the test results. However, our method uses hybrid degradation instead of linear degradation, and the load forces drop sharply after reaching the maximum load.

6.3. Discussion

The comparison of results suggests that the three criteria are equally effective in predicting the onset of damage. However, after damage initiation, the driving force of damage propagation differs for each criterion. For Max Stress and Shokreih-Hashin criteria, the axial compressive stress is the main factor, while for the proposed pw-FRP fiber kinking criterion, the driving force is the shear stress on the misaligned fiber frame. The proposed criterion is more compatible with the energy-based continuum stiffness degradation method because it accounts for fiber kinking damage by calculating the matrix damage around the misaligned fiber, and matrix damage propagation can be well-described by energy-based theories.

Using the proposed criteria with hybrid degradation, the damage evolution process is significantly shortened compared with linear degradation, even with the same fracture energy. The initial linear weakening stage of hybrid degradation allows for steady damage development and gradual structural stiffness loss, while the sudden failure at the end can well capture the abrupt structural collapse.

7. Conclusions

This paper conducted ultimate compression tests on two kinds of L-joints and examined the test phenomena, strain distribution, and load-displacement curve features. The test phenomenon revealed that the upper surface of the middle arc-shaped stiffener was the first to suffer damage under compression in the L-joints, and then the damage progressively expanded until the structure failed; the test data indicated that the structure could only maintain a linear elastic response under less than half of the failure load, and then the structure started to experience continuous stiffness degradation, but there was no evident damage phenomenon at this stage.

This paper then introduced a fiber kinking failure criterion and a hybrid stiffness reduction model for plain woven FRPs in simulation.

The proposed criterion determines whether fiber buckling failure occurs by checking the matrix damage state at the fiber bending location. For predicting damage initiation, the criterion can produce similar results as the existing conventional criteria. Moreover, its advantage lies in its integration with the energy based continuous damage evolution model. The fiber kinking criterion replaces the driving force of damage evolution from axial compressive stress to shear stress. As a result, the heavy fluctuations of the simulated L-D curves disappear, and the damage evolution process becomes much smoother as in the test.

The hybrid stiffness reduction model is adopted to simulate the damage evolution process. This model introduces a compulsory failure mechanism when the damage variable reaches a certain threshold. It overcomes the issue of slow damage propagation in the conventional continuous model without changing the fracture energy. It should be noted that the value of the threshold used in this paper is based on the author’s experience, which may require further investigation through experiments or mesoscopic simulation.

The proposed fiber kinking criterion and the hybrid degradation model are combined to predict the failure behavior of L-joints with different transitional radii. The simulation results show good agreement with the experimental data and reveal some features of the failure process that other theories fail to capture. The results showed that this method was more accurate and reliable than other methods, as it captured the failure behavior of both L-joints consistently.