Failure of blunt notched fiber metal laminates

Products: Abaqus/Standard  Abaqus/Explicit  

Fiber metal laminates (FMLs) are composed of laminated thin aluminum layers bonded with intermediate glass fiber-reinforced epoxy layers. FMLs are of great interest in the aerospace industry due to their superior properties, such as high fracture toughness and low-density when compared to solid aluminum sheets.

This example simulates failure and damage in a FML containing a blunt notch subjected to quasi-static loading conditions. Cohesive elements are used to model the interlaminar delamination, and the Abaqus damage model for fiber-reinforced materials is used to predict behavior of the fiber-reinforced epoxy layer. In addition, the behavior of the fiber-reinforced epoxy layer is also described using the model proposed by Linde et al. (2004), which is implemented in user subroutine UMAT. Both Abaqus/Standard and Abaqus/Explicit are used for simulation when the Abaqus built-in damage model is used for fiber-reinforced epoxy layers. This type of problem is important in the aerospace industry since blunt notches (e.g., fastener holes) commonly occur in airplane structures; the strength of the structure containing a blunt notch is a crucial design parameter. The models presented in this example demonstrate how to predict the blunt notch strength, the failure patterns of the fiber and matrix within the fiber-reinforced epoxy layer, and the delamination between different layers of FMLs.

Problem description and material characteristics

Figure 1.4.6�1 shows the geometry of the laminate containing the blunt notch for this example. The laminate is subjected to uniaxial tension in the longitudinal direction. The laminate is made of three layers of aluminum and two layers of 0�/90� glass fiber-reinforced epoxy. Only 1/8 of the laminate needs to be modeled, with appropriate symmetric boundary conditions applied as shown in Figure 1.4.6�2. Figure 1.4.6�2 also shows the through-thickness lay-up of the 1/8 model.

The material behavior of aluminum is assumed to be isotropic elastic-plastic with isotropic hardening. The Young’s modulus is 73800 MPa, and the Poisson’s ratio is 0.33; the isotropic hardening data are listed in Table 1.4.6�1.

The material behavior of the glass fiber-reinforced epoxy layers is assumed to be orthotropic, with stiffer response along the fiber direction and softer behavior in the matrix. The elastic properties—longitudinal modulus, ; transverse modulus, ; shear moduli,  and ; and Poisson’s ratios,  and —are listed in Table 1.4.6�2. The subscript “L” refers to the longitudinal direction (or fiber direction), and the subscript “T” refers to the two transverse directions orthogonal to the fiber direction. The damage initiation and evolution behavior is also assumed to be orthotropic. Table 1.4.6�3 lists the ultimate values of the longitudinal failure stresses,  and ; transverse failure stresses,  and ; and in-plane shear failure stress, . The superscripts “t” and “c” refer to tension and compression, respectively. The fracture energies of the fiber and matrix are assumed to be =12.5 N/mm and =1.0 N/mm, respectively.

Two material models that use the parameters listed above are considered, as follows:

1. The material is modeled based on the built-in model for damage in fiber-reinforced composites available in Abaqus (see “Damage and failure for fiber-reinforced composites: overview,” Section 23.3.1 of the Abaqus Analysis User's Manual).
2. The material is modeled using an alternative damage model that is based on the model proposed by Linde et al. (2004). The alternative damage model is implemented in user subroutine UMAT and is referred to in this discussion as the UMAT model. Details of the UMAT model are provided below.

The adhesive used to bond neighboring layers is modeled using interface layers with a thickness of t=0.001 mm. To simulate the interlaminar delamination, these interface layers are modeled with cohesive elements. The initial elastic properties of each interface are assumed to be isotropic with Young’s modulus E=2000 MPa and Poisson’s ratio =0.33. The failure stresses of the interface layers are assumed to be ===50 MPa; the fracture energies are ===4.0 N/mm. The subscripts “n,” “s,” and “t” refer to the normal direction and the first and second shear directions (for further discussion of the constitutive modeling methods used for the adhesive layers, see “Defining the constitutive response of cohesive elements using a traction-separation description,” Section 31.5.6 of the Abaqus Analysis User's Manual).

The plate is loaded with displacement boundary conditions applied at the right edge. To simplify the postprocessing, the displacement loading is applied at a reference point and an equation constraint is used to constrain the displacement along the loading direction between the right edge and the reference point. Except for those files designed exclusively to study the effect of the loading direction on the strength, the loading direction (along the global X-direction) aligns with the fiber direction of the 0� fiber-reinforced epoxy layer.

UMAT model for fiber-reinforced epoxy layers

For fiber-reinforced epoxy layers, the primary model considered is based on the built-in damage model for fiber-reinforced composites available in both Abaqus/Standard and Abaqus/Explicit. Alternatively, in Abaqus/Standard, the damage in the fiber-reinforced epoxy is also simulated using the model proposed by Linde et al. (2004), which is implemented in user subroutine UMAT and is discussed below.

In the UMAT model, the damage initiation criteria are expressed in terms of strains. Unlike the built-in model in Abaqus, which uses four internal (damage) variables, the UMATmodel uses two damage variables to describe damage in the fiber and matrix without distinguishing between tension and compression. Although the performance of the two models is expected to be similar for monotonic loads, such as in this example problem, the results obtained might differ considerably for more complex loads in which, for example, tension is followed by compression. For the UMAT model, if the material is subjected to tensile stresses that are large enough to cause partial or full damage (the damage variable corresponding to this damage mode will be greater than zero), both tensile and compressive responses of the material will be affected. However, in the case of the built-in damage model, only the tensile response will be degraded while the material compressive response will not be affected. In many cases the latter behavior is more suitable for modeling fiber-reinforced composites. In this section the governing equations for damage initiation and evolution as proposed by Linde et al. (2004) are discussed, followed by a description of the user subroutine UMAT implementation.

Damage in the fiber is initiated when the following criterion is reached:

where , , and  are the components of the elasticity matrix in the undamaged state. Once the above criterion is satisfied, the fiber damage variable, , evolves according to the equation

where  is the characteristic length associated with the material point. Similarly, damage initiation in the matrix is governed by the criterion

where , , and . The evolution law of the matrix damage variable, , is

During progressive damage the effective elasticity matrix is reduced by the two damage variables  and , as follows:

The use of the fracture energy-based damage evolution law and the introduction of the characteristic length  in the damage evolution law help to minimize the mesh sensitivity of the numerical results, which is a common problem of constitutive models with strain softening response. However, since the characteristic length calculation is based only on the element geometry without taking into account the real cracking direction, some level of mesh sensitivity remains. Therefore, elements with an aspect ratio close to one are recommended (for a discussion of mesh sensitivity, see “Concrete damaged plasticity,” Section 22.6.3 of the Abaqus Analysis User's Manual).

In user subroutine UMAT the stresses are updated according to the following equation:

The Jacobian matrix can be obtained by differentiating the above equation:

The above Jacobian matrix is not symmetric; therefore, the unsymmetric equation solution technique is recommended if the convergence rate is slow.

To improve convergence, a technique based on viscous regularization (a generalization of the Duvaut-Lions regularization) of the damage variables is implemented in the user subroutine. In this technique we do not use the damage variables calculated from the aforementioned damage evolution equations directly; instead, the damage variables are “regularized” via the following equations:

where  and  are the matrix and fiber damage variables calculated according to the damage evolution laws presented above,  and  are the “regularized” damage variables used in the real calculations of the damaged elasticity matrix and the Jacobian matrix, and  is the viscosity parameter controlling the rate at which the regularized damage variables  and  approach the true damage variables  and .

To update the “regularized” damage variables at time , the above equations are discretized in time as follows:

From the above expressions it can be seen that

Therefore, the Jacobian matrix can be further formulated as follows:

Care must be exercised to choose an appropriate value for  since a large value of viscosity might cause a noticeable delay in the degradation of the stiffness. To estimate the effect of viscous regularization, the approximate amount of energy associated with viscous regularization is integrated incrementally in user subroutine UMAT by updating the variable SCD as follows:

where  is the damaged elasticity matrix calculated using the damage variables,  and ; and  is the damaged elasticity matrix calculated using the regularized damage variables,  and . To avoid unrealistic results due to viscous regularization, the above calculated energy (available as output variable ALLCD) should be small compared to the other real energies in the system, such as the strain energy ALLSE.

This user subroutine can be used with either three-dimensional solid elements or elements with plane stress formulations. In the user subroutine the fiber direction is assumed to be along the local 1 material direction. Therefore, when solid elements are used or when shell elements are used and the fiber direction does not align with the global X-direction, a local material orientation should be specified. The damage variables—, , , and —are stored as solution-dependent variables, which can be viewed in the Visualization module of Abaqus/CAE.

Finite element model

The finite element model uses a separate mesh for each of the respective layers shown in Figure 1.4.6�2: two aluminium layers, two fiber-reinforced epoxy layers, and three adhesive layers. While not required, a similar finite element discretization in the plane of the laminate, such as that shown in Figure 1.4.6�3, can be used for all layers.

Modeling considerations for aluminium layers

Due to the interactions with the fiber-reinforced epoxy layers, the stress state within the aluminum layers (especially surrounding the notch tip) cannot be approximated using the plane stress assumption. To model this three-dimensional plasticity stress state accurately, solid elements must be used for the aluminum layers. In Abaqus/Standard incompatible mode elements (C3D8I) are used since local bending might exist in the post-failure region surrounding the notch. For the Abaqus/Explicit analysis, reduced-integration elements (C3D8R) are used for modeling the aluminum layers.

Modeling considerations for glass fiber-reinforced epoxy layers

The plane stress assumption can be used safely within the fiber-reinforced epoxy layers; therefore, either solid elements or shell elements can be adopted for these layers. However, it is important to have an accurate representation of the through-thickness geometry to model the interface between the adhesive and the fiber-reinforced epoxy realistically. This is achieved most conveniently with solid elements or continuum shell elements instead of conventional shell elements. The damage model for fiber-reinforced materials is available only for elements with a plane stress formulation. Therefore, continuum shell elements are used with this model. Models are also included in which continuum elements (C3D8R or C3D8) are used along with user subroutine UMAT to model the fiber-reinforced epoxy layers.

Modeling considerations for adhesive layers

Cohesive elements (COH3D8) are used for the interface layers. The elastic response is defined in terms of a traction-separation law with uncoupled behavior between the normal and shear components. For convenience, a constitutive thickness of 1.0 mm is used so that we do not need to distinguish between the separation displacement and the nominal strain (NE). However, since the actual thickness is 0.001 mm, the diagonal terms in the elasticity matrix need to be scaled by the inverse of the actual thickness as follows:

The quadratic nominal strain criterion is used for the damage initiation:

The damage evolution is based on fracture energy with the quadratic power law for the mixed mode behavior and exponential softening behavior (see “Defining the constitutive response of cohesive elements using a traction-separation description,” Section 31.5.6 of the Abaqus Analysis User's Manual).

Results and discussion

Results for each analysis are discussed in the following sections.

Abaqus/Standard results

Damage to the fiber-reinforced epoxy plays a key role in the response for the loading considered. Figure 1.4.6�4 shows the load-displacement curve for the 0� loading direction for both of the damage models considered for the fiber-reinforced epoxy. The response shows a “bilinear” shape before the sudden loss of loading capacity; i.e., an initial linear curve representing the initial elastic region, a smoothly deflecting nonlinear curve representing the local plasticity, and a second linear curve representing the net section yielding. The effect of the element type was studied using the UMAT model and C3D8R, C3D8, and SC8R elements; and the results are summarized in Figure 1.4.6�5 and Table 1.4.6�5. The numerical results obtained using different element types and different damage models are similar and show a good agreement with the experimental results of De Vries (2001).

The fiber and matrix damage patterns in the 0� fiber-reinforced epoxy layer at the failure load are shown in Figure 1.4.6�6 and Figure 1.4.6�7 for the built-in damage model for fiber-reinforced materials and in Figure 1.4.6�8 and Figure 1.4.6�9 for the UMAT model. It can be seen that the fiber damage in the 0� fiber-reinforced epoxy layer propagates along the ligament above the blunt notch tip (i.e., orthogonal to the loading direction). Figure 1.4.6�10 shows the matrix damage in the 90� layer for the damage model of Linde et al. (2004). There is no fiber damage in the 90� fiber-reinforced epoxy layer prior to the sudden fracture. Interlaminar damage is most severe between the 0� fiber-reinforced epoxy layer and the aluminum layer. These observations are in agreement with the experimental results of De Vries (2001).

Figure 1.4.6�11 and Table 1.4.6�4 give the load-displacement results for different values of the viscosity parameter, , obtained using the built-in damage model for fiber-reinforced materials. The same results for the UMAT model are given in Figure 1.4.6�12 and Table 1.4.6�6. The smaller the viscosity, the more abrupt the failure and the smaller the failure strength. Although a viscosity of 0.001 seems to overestimate the failure strength by a few percent (Table 1.4.6�4 and Table 1.4.6�6), the convergence is noticeably improved; thus, a viscosity of 0.001 is used for all the other studies in this example. For the built-in damage model for fiber-reinforced materials, only the viscosity in the fiber direction was varied while the viscosity in the matrix direction was kept constant at 0.005. This improved convergence and did not markedly affect the results.

The effect of the loading direction on the blunt notch strength is studied using the three-dimensional element, C3D8R, with the UMAT model. Three tests are performed in which the local material orientations in the 0�/90� fiber-reinforced epoxy are rotated by an angle of 15�, 30�, and 45�, respectively. For example, for a loading angle of 15� the fiber orientation in the 0� fiber-reinforced epoxy layer would be at a 15� angle with respect to the X-direction, while the fiber orientation in the 90� fiber-reinforced epoxy layer would be at an angle of �75� with respect to the X-direction (Figure 1.4.6�13). As can be seen in Figure 1.4.6�14, strain hardening is smaller for the larger loading angles. As can be seen in Figure 1.4.6�15, the failure strength decreases with the increasing loading angle and reaches the minimum at the 45� loading angle (the response for even larger loading angles is expected to be approximately symmetric with respect to the 45� angle due to the symmetric nature of the 0�/90� fiber-reinforced epoxy layer). As stated by De Vries (2001), this is expected and reflects the poor shear properties of the fiber-reinforced epoxy layer.

In the above discussions the net blunt notch strength is defined as , where  is the length of the ligament above the notch and t is the total thickness of the laminate. This example demonstrates that the approach employed in the study can be used to predict the blunt notch strength of the fiber metal laminates.

Abaqus/Explicit results

In the Abaqus/Explicit simulation we only consider loading along the 0� ply. The simulation is conducted without damage stabilization, and no mass scaling is used. However, in order to reduce the computational time, the total loading is applied in a short interval of time (0.001 s). The overall load-displacement curve obtained from the explicit dynamic simulation is compared with the Abaqus/Standard result (with viscosity of 0.001) in Figure 1.4.6�16. The results from the explicit dynamic simulation are presented using an antialiasing filter to remove high frequency noise (see “Filtering output and operating on output in Abaqus/Explicit” in “Output to the output database,” Section 4.1.3 of the Abaqus Analysis User's Manual). The overall response compares well with the Abaqus/Standard results with some differences in the peak value of the load and in the post-peak response. Note that damage stabilization is used in the Abaqus/Standard simulation to achieve convergence and is likely to change the overall response (especially in the post-peak portion of the load-displacement curve). On the other hand, the Abaqus/Explicit simulation does not use damage stabilization and is better able to capture the dynamic behavior inherent in the damage and failure processes. The contour plots of various damage variables in the 0� and 90� plies agree qualitatively with the corresponding plots obtained from the Abaqus/Standard simulation using the built-in damage model.

Python scripts

fml_c3d8r_deg0_vis1_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001.

fml_c3d8r_deg0_vis2_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.0004.

fml_c3d8r_deg0_vis3_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.00016.

fml_c3d8r_deg0_vis4_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.000064.

fml_c3d8_deg0_vis1_std.py

C3D8 used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001.

fml_sc8r_deg0_vis1_std.py

SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001.

fml_c3d8r_deg15_vis1_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 15�, and a viscosity of 0.001.

fml_c3d8r_deg30_vis1_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 30�, and a viscosity of 0.001.

fml_c3d8r_deg45_vis1_std.py

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 45�, and a viscosity of 0.001.

Input files

Abaqus/Standard input files

fml_frm_sc8r_deg0_vis001_std.inp

SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001 in the fiber direction (using built-in fiber-reinforced material damage model).

fml_frm_sc8r_deg0_vis0005_std.inp

SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.0005 in the fiber direction (using built-in fiber-reinforced material damage model).

fml_frm_sc8r_deg0_vis00025_std.inp

SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.00025 in the fiber direction (using built-in fiber-reinforced material damage model).

fml_c3d8r_deg0_vis1_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001 (using UMAT model).

fml_c3d8r_deg0_vis2_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.0004 (using UMAT model).

fml_c3d8r_deg0_vis3_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.00016 (using UMAT model).

fml_c3d8r_deg0_vis4_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.000064 (using UMAT model).

fml_c3d8_deg0_vis1_std.inp

C3D8 used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001 (using UMAT model).

fml_sc8r_deg0_vis1_std.inp

SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0�, and a viscosity of 0.001 (using UMAT model).

fml_c3d8r_deg15_vis1_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 15�, and a viscosity of 0.001 (using UMAT model).

fml_c3d8r_deg30_vis1_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 30�, and a viscosity of 0.001 (using UMAT model).

fml_c3d8r_deg45_vis1_std.inp

C3D8R used in the fiber-reinforced epoxy layer, a loading angle of 45�, and a viscosity of 0.001 (using UMAT model).

exa_fml_ortho_damage_umat.f

User subroutine UMAT for modeling the damage initiation and evolution in the fiber-reinforced epoxy layers.

Abaqus/Explicit input file

fml_frm_sc8r_deg0_exp.inp

SC8R elements used in the fiber-reinforced epoxy layer and a loading angle of 0� (using built-in fiber-reinforced material damage model).

References

De Vries,  T. J., “Blunt and Sharp Notch Behavior of Glare Laminates,” Ph.D dissertation, Delft University Press, 2001.

Hagenbeek,  M., C. Van Hengel, O. J. Bosker, and C. A. J. R. Vermeeren, “Static Properties of Fibre Metal Laminates,” Applied Composite Materials, vol. 10 207–222, 2003.

Linde,  P., J. Pleitner, H. De Boer, and C. Carmone, “Modelling and Simulation of Fiber Metal Laminates,” ABAQUS Users’ Conference, 2004.

Tables

Table 1.4.6�1 Isotropic hardening data for aluminum.

Yield stress (MPa)	Plastic strain (%)
300	0.000
320	0.016
340	0.047
355	0.119
375	0.449
390	1.036
410	2.130
430	3.439
450	5.133
470	8.000
484	14.710

Table 1.4.6�2 Orthotropic elastic properties of fiber-reinforced epoxy.

(MPa)	(MPa)	(MPa)	(MPa)
55000	9500	5500	3000	0.45	0.33

Table 1.4.6�3 Orthotropic damage initiation properties of fiber-reinforced epoxy.

(MPa)	(MPa)	(MPa)	(MPa)	(MPa)
2500	2000	50	150	50

Table 1.4.6�4 Net blunt notch strength (MPa) for different values of the viscosity parameter in fiber direction (using built-in fiber-reinforced material damage model, viscosity in the matrix direction =0.005).

Numerical results (SC8R, 0o loading angle)			Experimental results (De Vries, 2001)
=0.001	=0.0005	=0.00025
462.1	456.4	453.2	446

Table 1.4.6�5 Net blunt notch strength (MPa) for different element types used in the fiber-reinforced epoxy layers (using UMAT model).

Numerical results (=0.001, 0o loading angle)			Experimental results (De Vries, 2001)
C3D8R	C3D8	SC8R
463.7	467.1	458.7	446

Table 1.4.6�6 Net blunt notch strength (MPa) for different values of the viscosity parameter (using UMAT model).

Numerical results (C3D8R, 0o loading angle)				Experimental results (De Vries, 2001)
=0.001	=0.0004	=0.00016	=0.000064
463.7	453.8	449.2	448.2	446

Figures

Figure 1.4.6�1 Plate geometry.

Figure 1.4.6�2 (a) In-plane view of the 1/8 plate; (b) through-thickness lay-up of the 1/8 plate.

Figure 1.4.6�3 Finite element mesh.

Figure 1.4.6�4 Load-displacement curves for different damage models in the fiber-reinforced epoxy layer for the 0� loading direction, =0.001.

Figure 1.4.6�5 Load-displacement curves for different element types in the fiber-reinforced epoxy layer for the 0� loading direction (using UMAT model).

Figure 1.4.6�6 Fiber damage pattern in the 0� fiber-reinforced epoxy layer for the 0� loading direction (using built-in fiber-reinforced material model, DAMAGEFT contour plot).

Figure 1.4.6�7 Matrix damage pattern in the 0� fiber-reinforced epoxy layer for the 0� loading direction (using built-in fiber-reinforced material model, DAMAGEMT contour plot).

Figure 1.4.6�8 Fiber damage pattern in the 0� fiber-reinforced epoxy layer for the 0� loading direction (using UMAT model, SDV3 contour plot).

Figure 1.4.6�9 Matrix damage pattern in the 0� fiber-reinforced epoxy layer for the 0� loading direction (using UMAT model, SDV4 contour plot).

Figure 1.4.6�10 Matrix damage pattern in the 90� fiber-reinforced epoxy layer for the 0� loading direction (using UMAT model, SDV4 contour plot).

Figure 1.4.6�11 Load-displacement curves for different values of the viscosity parameter for the 0� loading direction (using built-in fiber-reinforced material damage model).

Figure 1.4.6�12 Load-displacement curves for different values of the viscosity parameter for the 0� loading direction (using UMAT).

Figure 1.4.6�13 Local material orientations in the fiber-reinforced epoxy layers for the 15� loading direction.

Figure 1.4.6�14 Load-displacement curves for different loading directions (using UMAT model).

Figure 1.4.6�15 Calculated blunt notch strength for different loading angles in comparison with the experimental results (using UMAT model).

Figure 1.4.6�16 Load-displacement curves for the 0� loading direction: Abaqus/Explicit versus Abaqus/Standard.