Influence of Geometric Nonlinearities on the Mixed Mode Decomposition in Asymmetric Dcb Samples

1. Introduction

Because fibre reinforced thermoset matrix composites are made up of sheets, delamination failure is one of the factors that limit the service life of these materials. Crack initiation and growth can take place in pure modes I, II and III. However, in most real situations a combination of these modes is usually observed.

The double cantilever beam (DCB) test configuration is usually used to obtain pure mode I and the end notch flexure (ENF) is commonly used to obtain pure mode II. These test procedures are elevated to the status of international standards [1,2,3].

On the other hand, there are several test procedures developed in the scientific literature to obtain pure mode III, although there is no scientific consensus to elevate any of them to international standard status. Some of these methods are the Split Cantilever Beam (SCB), Crack Rail Shear (CRS), Edge Crack Torsion test (ECT) and Longitudinal Semi Fixed Beam (LHFB) [4,5,6,7,8,9].

One of the most studied mixed modes in the scientific literature is mode I/II. The standardized Mixed Mode Bending (MMB) test method allows generating a mixed mode I/II at the crack tip [10,11]. In this test configuration, mode I is introduced by means of an opening movement, as in the DCB test, and mode II is produced by means of the specimen flexure (as in ENF test) (see Figure 1).

Another alternative test method to generate a mixed mode I/II is the Asymmetric Double Cantilever Beam test (ADCB) [12,13,14]. In this test configuration both sublaminates have different stiffness. This is the situation when we have sublaminates with different thickness and/or different elastic properties (e.g. bimaterial bonded joints). In this situation, a mixed mode I/II is generated at the crack tip. ADCB test specimens and fittings are as simple as in pure mode I tests (DCB). The ADCB test configuration allows generating mixed modes at the crack tip with predominance of mode I

Other test methods to obtain mixed modes I/II are based on the ENF test, as the Single Leg Bending (SLB) and Asymmetric Single Leg Bending (ASLB) tests [15,16].

There are other test configurations developed in the scientific literature to produce mixed mode I/II at the crack tip. This is the case of the Prestressed End-Notched Flexure (PENF) test (Figure 2) proposed by Szekrényes [17]. In this test configuration, mode I is introduced by means of a rod inserted between both sublaminates, and mode II is produced by the specimen flexure.

This test is simple to perform, as bonding hinges in the sample is not required. Boyano et al. studied the influence of different rod positions and diameters on the mixed mode at the crack tip [18].

The analytical determination of G_I and G_II is not simple when mixed mode is present at the crack tip. There are several analytical methods developed in the scientific literature to evaluate mode I/II partition, but there is not always consensus among researchers. [14,19,20,21,22].

The numerical calculation of mode I/II partition is usually carried out using the Virtual Crack Closure Technique (VCCT) [23] and the Cohesive Zone Modelling (CZM) [24,25] available in commercial FEM programs.

VCCT numerical calculations are usually carried out in global coordinates without taking into account the possible rotation of the specimen during the test. The same situation usually occurs with the analytical developments.

This article aims to draw attention to the fact that, if the ADCB specimen rotates in a non-negligible way during the test, the fact of not taking said rotation into account in the projection of the forces and displacements, can give rise to an appreciable error in the decomposition of modes I and II.

2. Materials and Methods

In this work, the ANSYS^® 2024 Academic Research package has been used to study by FEM the influence of the specimen rotation on the ERR decomposition in modes I and II (G_I and G_II) in ADCB specimens.

As mentioned above, the concept of modes I and II refers to the components of the energy release rate calculated in the direction perpendicular and tangential to the delamination plane respectively. To correctly calculate G_I and G_II, the components of the forces and displacements involved in the energy calculation must be projected in directions perpendicular and tangential to the delamination plane (x_L-y_L coordinate system) (see Figure 3).

In this study, two coordinate systems are defined: a global coordinate system (GCS) (x_G-y_G), with the x_G-axis aligned with the crack plane when the specimen is unloaded, that remain horizontal during the test, and a local coordinate system (LCS) (x_L-y_L) that rotates with the crack plane while loading the sample. Therefore, the x_L-axis is always aligned with the crack during the test (see Figure 3).

If the rotation of the specimen during the test is large, the difference in the calculations obtained using both coordinate systems may be non-negligible. As shown in Figure 3, during the ADCB test, when the sublaminates have a significant difference in stiffness, an appreciable rotation of the specimen occurs.

To carry out this study, several ADCB models were prepared with different thickness ratios (h₁/h₂) in order to obtain different degrees of modes I/II (see Figure 4).

A constant displacement δ was applied to all models on the sample lips. G_I and G_II were then calculated using forces and displacements in global and local coordinates.

FEM packages usually implement the Virtual Crack Closure Technique (VCCT) in order to calculate G_I and G_II. The VCCT method calculates the energy released by the opening of the nodes during crack growth in a single step. This is a simplification since the force is calculated at the crack tip when the crack has already grown, but it is a very efficient calculation method.

In the VCCT method [26] (see Figure 5), the components of forces and displacements are calculated in the same step using two different node pairs. Forces are calculated in nodes 2i-2i’ and displacements are calculated in nodes 1i-1i’. This approximation is valid providing that forces at the crack tip before and after the crack extension are similar when the crack increment (element size) is small enough.

The calculation using the VCCT method was performed with the CINT function of the Ansys program (CINT,TYPE,VCCT) [27].

G_I and G_II are then calculated as follows:

$$G_I={\displaystyle \frac{1}{2B∆a}}\sum _{i=1}^n{F}_{y2i}(v_{1i}-v_{1^{\text{'}}i})$$

(1)

$$G_{II}={\displaystyle \frac{1}{2B∆a}}\sum _{i=1}^n{F}_{x2i}\left(u_{1i}-u_{1^{\text{'}}i}\right)$$

(2)

B is the specimen width. The subscript i takes into account a 3D system with n nodes along the crack front.

The VCCT method, implemented in this commercial FEM package, uses the global coordinate system to perform de calculations (the coordinate system is not updated as the crack plane rotates).

In this article, a similar method without simplification (TSEM) is also used. In this method, loads are calculated in the first step and the corresponding displacements are calculated in a second step once the nodes have separated. Loads and displacements are therefore calculated in the same node pair (nodes 1i-1i’ in Figure 6 and Figure 7) [28,29].

G_I and G_II are calculated by means of TSEM as follows:

$$G_I={\displaystyle \frac{1}{2B∆a}}\sum _{i=1}^n{F}_{y1i}(v_{1i}-v_{1^{\text{'}}i})$$

(3)

$$G_{II}={\displaystyle \frac{1}{2B∆a}}\sum _{i=1}^n{F}_{x1i}\left(u_{1i}-u_{1^{\text{'}}i}\right)$$

(4)

Therefore, data needs to be collected in two successive calculation steps to perform the energy calculations. This method was applied using both global (Figure 6) and local coordinate systems (Figure 7).

This method is not originally implemented in the ANSYS FEM package. Therefore, a script programmed in APDL (Ansys Parametric Design Language) was specifically developed for this purpose in this work. This method was programmed using local and global coordinates. In this way, it will be possible to compare results taking or not into account the rotation of the coordinate system. Moreover, TSEM results in global coordinates can be compared with the VCCT originally implemented in Ansys (and most commercial FEM packages), as they are equivalent when the element size ahead of the crack tip is small enough.

The material used to perform this study is a unidirectional carbon/epoxy laminate AS4/8552. The mechanical properties of this laminate taken from [18] are shown in Table 1. Poisson coefficient ν_yz has been calculated as isotropic.

3. Model Calibration

A 2D FEM model was prepared using ANSYS^® 2024 Academic Research. The element used was the 2D structural PLANE182 with the plane strain option.

A first model was prepared to optimize the element size near the crack tip. This model had the following dimensions:

• L=150 mm
• a₀=50 mm
• h₁=4.6 mm
• h₂=1.4 mm

A regular mesh was implemented with the same size before and after the crack tip (Figure 8). Several runs using the TSEM procedure were performed in global coordinates with different element sizes to optimize the model. The mesh was refined around the crack tip until results variation was negligible. The elements were nearly squared.

A displacement of δ=8 mm was applied to the sample lips.

Finally, the element size in the vicinity of the crack tip was set to 0.100 mm for the following studies. A thinner mesh would produce a negligible change in ERR values (0.1%).

Another issue to take in consideration is whether the calculations should be performed in a linear or non-linear way due to sample rotation.

To perform VCCT crack growth simulations, the Ansys manual specifies the assumption that the model undergoes a small deformation (or a small rotation), and so, the calculation is performed with the option NLGEOM, OFF [27]. This option causes the calculation to be performed by a linear procedure, without taking into account possible non-linear effects. There are three possible sources of non-linearities: material non-linearity, contacts and large displacements or deformations. This is an issue to take into account, as asymmetric specimens, may show non negligible rotations.

Large displacements or rotations can cause the stiffness of the system to change as the load is applied to the sample lips. If the change in stiffness of the system during the loading process is not negligible, it must be taken into account as a source of non-linearity in the calculations.

Two runs were performed to investigate the influence of the calculation procedure (linear or non-linear procedure).

The model used was the same as the one used for element size optimisation. In this case, a displacement of 4.5 mm was applied to the sample. This displacement is assumed to be close to the critical displacement of the material for this configuration. The obtained results are shown in Table 2.

The units of the applied load were set as N assuming a width of 1 mm (a 2D model has been used). In this Table, the total energy release rate (G_T) was calculated in two ways: as the sum of G_I and G_II previously calculated by means of the TSEM procedure, and, directly, by means of the overall loss of elastic energy divided by the area of the crack increment (one element length). In mathematical form:

$$G_T=-{\displaystyle \frac{1}{B}}{\displaystyle \frac{dU}{da}}$$

(5)

Both procedures give similar results for both calculation procedures (linear and non-linear), but, as can be seen in the Table above, linear calculations give identical results irrespective of the procedure used to calculate G_T (0.382 N/mm).

This model was also used to compare both FEM solution procedures with the analytical Corrected Beam Theory proposed by ISO 15024 [1] standard, equivalent to the Modified Beam Theory in ASTM D 5528 [2]. Although this formulation is designed to calculate G_I in symmetric samples, it can be also used to calculate G_T in asymmetric samples as it was demonstrated in [30].

$$G_T={\displaystyle \frac{3P\delta }{2B\left(a+∆\right)}}$$

(6)

where,

• P: applied load
• δ: applied displacement
• B: sample width
• a: crack length
• Δ: crack length correction

Δ takes into account that the beam is not perfectly built in. This parameter is calculated by plotting the cube-root of the compliance as a function of delamination length. The extrapolation of the linear fit through the data in the plot yields Δ as the x-intercept [1]. This was done by analysing various models with h₁=4.6 mm, h₂=1.4 mm and different initial crack lengths by FEM. Once Δ was obtained, a new model with a₀= 50 mm and an applied displacement δ=4.5 mm was analysed. Table 3 shows the obtained results calculated following equations 3 and 4, and additionally, with the analytical equation 6. These calculations were performed using linear and non-linear solutions.

As can be seen in this Table, both solution procedures (linear and non-linear) yield similar results and agree very well with the analytical solution. The time needed to solve the model with each procedure is very similar since a 2D model is used.

Because the solution time is similar between the two procedures, non-linear solution will be used in the following calculations to take into account the small non-linearities present in this model (the stiffness matrix of the model changes slightly as the displacement is applied).

4. Influence of the Sample Rotation

Once the base model has been optimized, several runs were performed increasing the applied displacement with the same model used in the previous studies. This causes the sample to progressively increase its rotation (see Figure 9).

These models were analysed using TSEM in local and global coordinates by means of the APDL script developed specifically for this purpose and also, using the VCCT procedure available in ANSYS by means of the CINT,TYPE,VCCT function. As commented before, this function, implemented in the Ansys package, uses only global coordinates. Results are shown in Table 4, Table 5 and Table 6.

Comparing Table 4 and Table 5, we can see that both procedures yield practically identical results and match the partition values G_I/G_T and G_II/G_T. These results were expected as both procedures perform the calculations in global coordinates, and they are very similar: only differs in the node where the force is calculated (one element length away from each other). The two procedures tend to converge on each other as the size of the element at the crack tip decreases.

We can also observe in Table 4 and Table 5, that as the rotation of the sample increases, the partition values (G_I/G_T - G_II/G_T) change from 75%-25% for δ=1 mm to 61%-39% for δ=20 mm. This represents a change of 19% in the G_I/G_T values and 58% in the G_II/G_T values. This change in the partition values can be attributed to the appreciable rotation of the sample as the displacement increases. We must remember that, in these tables, forces and displacements are calculated always horizontal and vertical (according to the GCS) independently of the sample rotation.

On the other hand, Table 6 shows the calculations using forces and displacements always perpendicular and parallel to the crack plane, this is, considering the rotation of the sample (LCS).

This procedure gives rise to different results from the previous ones using GCS and the difference between results increases as the rotation of the specimen increases. Using local coordinates, G_I/G_T - G_II/G_T ratios change from 76%-24% for δ=1 mm to 79%-21% for δ=20 mm as the displacement increases. As we can see, the ERR partition has a much smaller change as the sample increases rotation. In this case, the small change in the partition values could be influenced, among other issues, by the shortening of the distance a₀ due to sample rotation and arm flexure.

Finally, as it was expected, the total energy G_T dissipated when the crack grows is the same for a given displacement, independently of the coordinate system used (Table 5 and Table 6). Obviously, the total energy dos not dependent on the coordinate system used to calculate G_I and G_II.

We can evaluate ERR modes I and II in global coordinates (G^G_I, G^G_II) as a function of modes I and II in local coordinates (G^L_I, G^L_II) (Figure 10) or vice versa as follows:

$$F_y^G=F_y^L\cos \left(\alpha \right)+F_x^L\sin \left(\alpha \right), { u}_y^G=u_y^L\cos \left(\alpha \right)+u_x^L\sin \left(\alpha \right)$$

(7)

$$F_x^G=F_x^L\cos \left(\alpha \right)-F_y^L\sin \left(\alpha \right), u_x^G=u_x^L\cos \left(\alpha \right)-u_y^L\sin \left(\alpha \right)$$

(8)

$$G_I^G=G_I^L{cos}^2\alpha +G_{II}^L{sin}^2\alpha +{\displaystyle \frac{1}{2B∆a}}\left(F_y^L\mathsf{\Delta }{u}_x^L+F_x^L\mathsf{\Delta }{u}_y^L\right)\sin \left(\alpha \right)\cos \left(\alpha \right)$$

(9)

$$G_{II}^G=G_I^L{sin}^2\alpha +G_{II}^L{cos}^2\alpha -{\displaystyle \frac{1}{2B∆a}}\left(F_y^L\mathsf{\Delta }{u}_x^L+F_x^L\mathsf{\Delta }{u}_y^L\right)\sin \left(\alpha \right)\cos \left(\alpha \right)$$

(10)

and adding equations 9 and 10:

$$G_T=G_I^G+G_{II}^G=G_I^L+G_{II}^L$$

(11)

In Table 7 we can see the differences between the calculated values of G_I and G_II using GCS and LCS for different displacements applied to the sample (calculated based on results in Table 5 and Table 6).

As expected, we can see that, for small rotations the difference between calculations in GCS and LCS is low, but if the rotation of the sample is significant, the error made when calculating in GCS can be very significant compared to LCS. For a displacement of only 5 mm, the difference obtained in G_I and G_II values is 4% and 14% respectively when calculating in GCS and LCS. And for a displacement of 10 mm, de differences in G_I and G_II values reach 10% and 35% respectively. This should be considered in the numerical and analytical calculations when the sample rotation is large. That is, when the applied displacement, and/or the difference in the stiffness of both sublaminates is considerable.

As stated above, another source of non-linearity is the displacement of the crack tip node as the applied load increases (see Table 7). For δ=20 mm, the equivalent crack length (measured perpendicular to the loading axis) shortens 1.8 mm making the sample response stiffer.

Next, some models were prepared with different h₂/h₁ and E₂/E₁ ratios to produce different mixity-modes. The applied load was set to δ=10 mm and the initial crack length was set to a₀=50 mm. The obtained results can be seen in Table 8.

The ratio E₁h₁²/E₂h₂² has been included in this Table to give an idea of how far each configuration is from pure mode I condition, assuming that pure mode I condition is reached when the longitudinal strain distributions at both sublaminates interfaces are identical [31] (Figure 11 and Figure 12). Therefore, pure mode I condition is achieved as follows:

$$u_1=u_2\to {\gamma }_{xy}=0$$

(12)

$${\displaystyle \frac{d{u}_1}{dx}}={\displaystyle \frac{d{u}_2}{dx}}\to {\epsilon }_1={\epsilon }_2\to {\displaystyle \frac{M{h}_1}{2E_1{I}_1}}={\displaystyle \frac{M{h}_2}{2E_2{I}_2}}\to E_1{h}_1^2=E_2{h}_2^2$$

(13)

Figure 13 shows the rotation of the sample with h₂/h₁=0.2 compared to the symmetric sample for δ=20 mm.

Finally, crack growth calculations were performed in order to evaluate the influence of the method used (LCS or GCS) on the obtained critical load.

To determine the crack onset and subsequent crack growth, it is necessary to define a fracture criterion. The Benzeggagh-Kenane criterion was assumed [27,32]:

$${\displaystyle \frac{G_T}{G_{Ic}+\left(G_{IIc}-G_{Ic}\right){\left(\frac{G_{II}}{G_T}\right)}^{\eta }}}\ge 1$$

(14)

The exponent η was set to 1.45 following [33]

The model was set as follows:

• E₁=E₂=144,000 MPa
• h₁= 5 mm
• h₂= 1 mm
• a₀= 50 mm

The obtained results are shown in Table 9. As can be seen in this Table, the difference between critical values (both displacement and load) is in the order of 4% in this case. The calculation performed at LCS is on the safe side as it predicts a crack growth onset at lower load and displacement. Figure 14 shows the corresponding load-displacement curves.

5. Conclusions

The ADCB configuration is, essentially, a DCB specimen with different sublaminate thickness and/or different elastic modulus. This is the case, for example, in bonded joints with dissimilar adherents.

When the ADCB specimen suffers large rotations during the application of the load, this nonlinearity must be considered in the calculation of displacements and forces, both in numerical and analytical procedures. In this work it has been demonstrated that, the difference in the ERR calculations when the angle of rotation of the sample is taken into account or not, can be very significant depending on the sample configuration and the applied displacement to the sample lips.

For a displacement of only 5 mm, the difference obtained in G_I and G_II values is 4% and 14% respectively when calculating in GCS and LCS. And for a displacement of 10 mm, de differences in G_I and G_II values reach 10% and 35% respectively.

When performing crack growth calculations, the difference between critical values (both displacement and load at the crack onset) can be in the order of 4%.

In the case of large rotation, it must be ensured that the ERR in modes I and II is calculated using the displacements and forces projected according to a coordinate system that rotates with the crack plane, that is, perpendicular and tangential to the crack plane. This is the case when the applied displacement is large, and/or when the difference in stiffness of the sub-laminates is large.

This correction is not usually implemented in numerical codes or in analytical developments, so this is an issue to take into account when the rotation of the sample is not negligible.