Failure of blunt notched fiber metal laminates

This example simulates failure and damage in a fiber metal laminate containing a blunt notch subjected to quasi-static loading conditions.

This page discusses:

Problem description and material characteristics
UMAT model for fiber-reinforced epoxy layers
Multiscale model for fiber-reinforced epoxy layers
Finite element model
Results and discussion
Python scripts
Input files
References
Tables
Figures

ProductsAbaqus/StandardAbaqus/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.

Cohesive elements are used to model the interlaminar delamination, and three different material models are used to predict the behavior of the fiber-reinforced layer:

Hashin damage model for unidirectional fiber-reinforced materials (Hashin Criterion).
Damage model proposed by Linde et al. (2004), which is implemented in user subroutine UMAT.
Abaqus multiscale material model (Multiscale Material Modeling).

Both Abaqus/Standard and Abaqus/Explicit are used for simulation when the Hashin damage model is used for the fiber-reinforced epoxy layers.

This type of problem is important in the aerospace industry since blunt notches (for example, 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 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 2. Figure 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 73,800 MPa, and the Poisson’s ratio is 0.33; the isotropic hardening data are listed in Table 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, $E_{L}$ ; transverse modulus, $E_{T}$ ; shear moduli, $G_{L T}$ and $G_{T T}$ ; and Poisson’s ratios, $ν_{T T}$ and $ν_{L T}$ —are listed in Table 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 3 lists the ultimate values of the longitudinal failure stresses, $σ_{L}^{f, t}$ and $σ_{L}^{f, c}$ ; transverse failure stresses, $σ_{T}^{f, t}$ and $σ_{T}^{f, c}$ ; and in-plane shear failure stress, $τ_{L T}^{f}$ . The superscripts “t” and “c” refer to tension and compression, respectively. The fracture energies of the fiber and matrix are assumed to be $G_{f}$ =12.5 N/mm and $G_{m}$ =1.0 N/mm, respectively.

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

The material is modeled based on the Hashin model for damage in unidirectional fiber-reinforced composites available in Abaqus (see About Damage and Failure for Fiber-Reinforced Composites).
The material is modeled using the damage model proposed by Linde et al. (2004). This 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 material is modeled using the multiscale material model in Abaqus (see Mean-Field Homogenization), and the damage model is implemented in user subroutine USDFLD at the constituent level. Details of the multiscale material and the damage 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 $t_{n}^{f}$ = $t_{s}^{f}$ = $t_{t}^{f}$ =50 MPa; the fracture energies are $G_{n}$ = $G_{s}$ = $G_{t}$ =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).

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 Hashin damage model for unidirectional 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 Hashin damage model in Abaqus, which uses four internal (damage) variables, the UMAT model 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 is greater than zero), both tensile and compressive responses of the material are affected. However, in the case of the Hashin damage model, only the tensile response is degraded while the material compressive response is not 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:

f_{f} = \sqrt{\frac{ϵ_{11}^{f, t}}{ϵ_{11}^{f, c}} {(ϵ_{11})}^{2} + (ϵ_{11}^{f, t} - \frac{{(ϵ_{11}^{f, t})}^{2}}{ϵ_{11}^{f, c}}) ϵ_{11}} > ϵ_{11}^{f, t},

where $ϵ_{11}^{f, t} = σ_{L}^{f, t} / C_{11}$ , $ϵ_{11}^{f, c} = σ_{L}^{f, c} / C_{11}$ , and $C_{i j}$ are the components of the elasticity matrix in the undamaged state. Once the above criterion is satisfied, the fiber damage variable, $d_{f}$ , evolves according to the equation

d_{f} = 1 - \frac{ϵ_{11}^{f, t}}{f_{f}} e^{(- C_{11} ϵ_{11}^{f, t} (f_{f} - ϵ_{11}^{f, t}) L^{c} / G_{f})},

where $L^{c}$ is the characteristic length associated with the material point. Similarly, damage initiation in the matrix is governed by the criterion

f_{m} = \sqrt{\frac{ϵ_{22}^{f, t}}{ϵ_{22}^{f, c}} {(ϵ_{22})}^{2} + (ϵ_{22}^{f, t} - \frac{{(ϵ_{22}^{f, t})}^{2}}{ϵ_{22}^{f, c}}) ϵ_{22} + {(\frac{ϵ_{22}^{f, t}}{ϵ_{12}^{f}})}^{2} {(ϵ_{12})}^{2}} > ϵ_{22}^{f, t},

where $ϵ_{22}^{f, t} = σ_{T}^{f, t} / C_{22}$ , $ϵ_{22}^{f, c} = σ_{T}^{f, c} / C_{22}$ , and $ϵ_{12}^{f} = τ_{L T}^{f} / C_{44}$ . The evolution law of the matrix damage variable, $d_{m}$ , is

d_{m} = 1 - \frac{ϵ_{22}^{f, t}}{f_{m}} e^{(- C_{22} ϵ_{22}^{f, t} (f_{m} - ϵ_{22}^{f, t}) L^{c} / G_{m})} .

During progressive damage the effective elasticity matrix is reduced by the two damage variables $d_{f}$ and $d_{m}$ , as follows:

C_{d} = [\begin{matrix} (1 - d_{f}) C_{11} & (1 - d_{f}) (1 - d_{m}) C_{12} & (1 - d_{f}) C_{13} & 0 & 0 & 0 \\ (1 - d_{m}) C_{22} & (1 - d_{m}) C_{23} & 0 & 0 & 0 \\ C_{33} & 0 & 0 & 0 \\ s y m m e t r i c & (1 - d_{f}) (1 - d_{m}) C_{44} & 0 & 0 \\ C_{55} & 0 \\ C_{66} \end{matrix}] .

The use of the fracture energy-based damage evolution law and the introduction of the characteristic length $L^{c}$ 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).

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

σ = C_{d} : ϵ .

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

\begin{matrix} \frac{\partial σ}{\partial ϵ} & = C_{d} + \frac{\partial C_{d}}{\partial ϵ} : ϵ \\ = C_{d} + (\frac{\partial C_{d}}{\partial d_{m}} : ϵ) (\frac{\partial d_{m}}{\partial f_{m}} \frac{\partial f_{m}}{\partial ϵ}) + (\frac{\partial C_{d}}{\partial d_{f}} : ϵ) (\frac{\partial d_{f}}{\partial f_{f}} \frac{\partial f_{f}}{\partial ϵ}) . \end{matrix}

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” using the following equations:

{\dot{d}}_{m}^{v} = \frac{1}{η} (d_{m} - d_{m}^{v}),

{\dot{d}}_{f}^{v} = \frac{1}{η} (d_{f} - d_{f}^{v}),

where $d_{m}$ and $d_{f}$ are the matrix and fiber damage variables calculated according to the damage evolution laws presented above, $d_{m}^{v}$ and $d_{f}^{v}$ 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 $d_{m}^{v}$ and $d_{f}^{v}$ approach the true damage variables $d_{m}$ and $d_{f}$ .

To update the “regularized” damage variables at time $t_{0} + Δ t$ , the above equations are discretized in time as follows:

{d_{m}^{v} |}_{t_{0} + Δ t} = {\frac{Δ t}{η + Δ t} d_{m} |}_{t_{0} + Δ t} + {\frac{η}{η + Δ t} d_{m}^{v} |}_{t_{0}},

{d_{f}^{v} |}_{t_{0} + Δ t} = {\frac{Δ t}{η + Δ t} d_{f} |}_{t_{0} + Δ t} + {\frac{η}{η + Δ t} d_{f}^{v} |}_{t_{0}} .

From the above expressions, it can be seen that

\frac{\partial d_{m}^{v}}{\partial d_{m}} = \frac{\partial d_{f}^{v}}{\partial d_{f}} = \frac{Δ t}{η + Δ t} .

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

\frac{\partial Δ σ}{\partial Δ ϵ} = C_{d} + [(\frac{\partial C_{d}}{\partial d_{m}^{v}} : ϵ) (\frac{\partial d_{m}}{\partial f_{m}} \frac{\partial f_{m}}{\partial ϵ}) + (\frac{\partial C_{d}}{\partial d_{f}^{v}} : ϵ) (\frac{\partial d_{f}}{\partial f_{f}} \frac{\partial f_{f}}{\partial ϵ})] \frac{Δ t}{η + Δ t} .

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:

Δ E_{S C D} = \frac{1}{2} [(C_{d} : ϵ) |_{t_{0}} + (C_{d} : ϵ) |_{t_{0} + Δ t}] : d ϵ - \frac{1}{2} [(C_{d}^{0} : ϵ) |_{t_{0}} + (C_{d}^{0} : ϵ) |_{t_{0} + Δ t}] : d ϵ,

where $C_{d}^{0}$ is the damaged elasticity matrix calculated using the damage variables, $d_{m}$ and $d_{f}$ ; and $C_{d}$ is the damaged elasticity matrix calculated using the regularized damage variables, $d_{m}^{v}$ and $d_{f}^{v}$ . 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— $d_{m}$ , $d_{f}$ , $d_{m}^{v}$ , and $d_{f}^{v}$ —are stored as solution-dependent variables.

Multiscale model for fiber-reinforced epoxy layers

The material of the fiber-reinforced epoxy layers is also modeled with the multiscale material model available in Abaqus/Standard. The Mori-Tanaka homogenization method is used, and the volume fraction of the fiber is set to 50%, which is typical for fiber-reinforced epoxy composites. Both the matrix material and fiber material are modeled with a linear elastic material model. The elastic moduli are calibrated to match the composite properties listed in Table 2. The matrix material is assumed to be isotropic. The Poisson's ratio is 0.41, and the calibrated Young's modulus is 6161.4 MPa. To better match the composite properties, the fiber material is modeled with an orthotropic linear elastic material with the calibrated properties listed in Table 4. The shape of the fiber is assumed to be cylindrical with infinite length.

In the multiscale material model, the damage initiation criteria are expressed at the constituent level rather than at the composite level. The stress limit values of the matrix material and fiber material are obtained from microlevel solutions when the stress limits listed in Table 3 are reached at the composite level.

Damage in the fiber material is initiated when one of the following criteria is reached:

Fiber tension $(σ_{11}^{f} > 0)$ :
$F_{f}^{t} = {(\frac{σ_{11}^{f}}{X_{f}^{T}})}^{2} \geq 1.$
Fiber compression $(σ_{11}^{f} < 0)$ :
$F_{f}^{t} = {(\frac{σ_{11}^{f}}{X_{f}^{C}})}^{2} \geq 1.$

Damage in the matrix material is initiated when one of the following criteria is reached:

Matrix tension $(σ_{22}^{m} > 0)$ :
$F_{m}^{t} = {(\frac{σ_{22}^{m}}{Y_{m}^{T}})}^{2} + {(\frac{σ_{12}^{m}}{S_{m}^{L}})}^{2} \geq 1.$
Matrix compression $(σ_{22}^{m} < 0)$ :
$F_{m}^{c} = {(\frac{σ_{22}^{m}}{Y_{m}^{C}})}^{2} + {(\frac{σ_{12}^{m}}{S_{m}^{L}})}^{2} + {(\frac{σ_{23}^{m}}{S_{m}^{T}})}^{2} \geq 1.$

The parameters in the above equations are listed in Table 5. The maximum values of the failure indices $F_{f}^{T}$ , $F_{f}^{C}$ , $F_{m}^{T}$ , and $F_{m}^{C}$ are saved as solution-dependent state variables. Once the maximum value of the failure index exceeds 1.0, the damage variable $d_{f, m}$ is set to 1.0 and continues to have the value of 1.0 even though the stresses might reduce significantly, which ensures that the material does not "heal" after it is damaged. To improve convergence, the damage variables are "regularized" using the following equations:

{\dot{d}}_{m}^{v} = \frac{1}{η_{m}} (d_{m} - d_{m}^{v}),

{\dot{d}}_{f}^{v} = \frac{1}{η_{f}} (d_{f} - d_{f}^{v}) .

The value of the "regularized" damage variable is assigned to the field variable used to define the elastic properties of the matrix and fiber materials. The Poisson's ratios remain unchanged while the other elastic moduli reduce by 1×10⁻⁶ when the field variables reach 1.0.

Both damage initiation and the viscous regularization of the damage variables are implemented in user subroutine USDFLD.

Finite element model

The finite element model uses a separate mesh for each of the respective layers shown in Figure 2: two aluminum 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 3, can be used for all layers.

Modeling considerations for aluminum 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 Hashin damage model for unidirectional 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. The multiscale material model using continuum shell elements appears to give the best results for the damage model described above. The model is also available for plane stress elements and three-dimensional continuum solid elements.

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:

K_{n n} = \frac{E}{t}, K_{s s} = \frac{G}{t}, K_{t t} = \frac{G}{t} .

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

\sqrt{{(\frac{ϵ_{n}}{ϵ_{n}^{f}})}^{2} + {(\frac{ϵ_{s}}{ϵ_{s}^{f}})}^{2} + {(\frac{ϵ_{t}}{ϵ_{t}^{f}})}^{2}} > 1.0, ϵ_{n}^{f} = \frac{t_{n}^{f}}{E / t}, ϵ_{s}^{f} = \frac{t_{s}^{f}}{G / t}, ϵ_{t}^{f} = \frac{t_{t}^{f}}{G / t} .

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).

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 4 shows the load-displacement curve for the 0° loading direction for the Linde and Hashin (built-in) the damage models considered for the fiber-reinforced epoxy. The response shows a “bilinear” shape before the sudden loss of loading capacity; that is, 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 5 and Table 7. 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 6 and Figure 7, for the Hashin damage model for fiber-reinforced materials;
Figure 8 and Figure 9, for the UMAT model; and
Figure 10 and Figure 11, for the multiscale material 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 (that is, orthogonal to the loading direction).

Figure 12 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).

The load-displacement results for different values of the viscosity parameter, $η$ are given in

Figure 13 and Table 6, for the Hashin damage model for fiber-reinforced materials;
Figure 14 and Table 8, for the UMAT model; and
Figure 15, for the multiscale material model.

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 6 and Table 8), the convergence is noticeably improved. Therefore, a viscosity of 0.001 is used for all the other studies in this example.

For the Hashin 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. For the multiscale material model, the viscosity in the matrix material was also kept at 0.005, while the viscosity in the fiber was varied between 0.0025 and 0.002. Since the viscosity parameter in this model also dictates the damage evolution, it has a larger impact on the failure strength, as shown in Table 9.

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 16). As can be seen in Figure 17, strain hardening is smaller for the larger loading angles. As can be seen in Figure 18, 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 $P_{m a x} / l_{1} t$ , where $l_{1}$ 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 19. 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). 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. 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 Hashin 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 Hashin 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 Hashin 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 Hashin 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.
fml_mfh_sc8r_deg0_vis002_std.inp: SC8R used in the fiber-reinforced epoxy layer, a loading angle of 0°, and a viscosity of 0.002 in the fiber material (using multiscale material model).
fml_mfh_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 material (using multiscale material model).
fml_mfh_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 material (using multiscale material model).
fml_mfh_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 material (using multiscale material model).
fml_material_mfh_dmgini.inp: Include file with material definition of the multiscale material.
exa_fml_mfh_damage_usdfld.f: User subroutine USDFLD for modeling the damage initiation and evolution in the fiber and matrix material (using multiscale material model).

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 Hashin 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. 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 2. Orthotropic elastic properties of fiber-reinforced epoxy.
$E_{L}$ (MPa)	$E_{T}$ (MPa)	$G_{L T}$ (MPa)	$G_{T T}$ (MPa)	$ν_{T T}$	$ν_{L T}$
55000	9500	5500	3000	0.45	0.33

Table 3. Orthotropic damage initiation properties of fiber-reinforced epoxy.
$σ_{L}^{f, t}$ (MPa)	$σ_{L}^{f, c}$ (MPa)	$σ_{T}^{f, t}$ (MPa)	$σ_{T}^{f, c}$ (MPa)	$τ_{L T}^{f}$ (MPa)
2500	2000	50	150	50

Table 4. Orthotropic elastic properties of fiber material.
$E_{L}$ (MPa)	$E_{T}$ (MPa)	$G_{L T}$ (MPa)	$G_{T T}$ (MPa)	$ν_{T T}$	$ν_{L T}$
103739	12170	29655.4	4261.56	0.242	0.1594

Table 5. Damage initiation properties of matrix and fiber in the multiscale material.
$X_{f}^{T}$ (MPa)	$X_{f}^{C}$ (MPa)	$Y_{m}^{T}$ (MPa)	$Y_{m}^{C}$ (MPa)	$S_{m}^{L}$ (MPa)	$S_{m}^{T}$ (MPa)
4550	3700	48	142	30	66

Table 6. Net blunt notch strength (MPa) for different values of the viscosity parameter in fiber direction (using Hashin damage model, viscosity in the matrix direction $η_{m}$ =0.005).
Numerical results (SC8R, 0^o loading angle)			Experimental results(De Vries, 2001)
$η_{f}$ =0.001	$η_{f}$ =0.0005	$η_{f}$ =0.00025	Experimental results(De Vries, 2001)
462.1	456.4	453.2	446

Table 7. Net blunt notch strength (MPa) for different element types used in the fiber-reinforced epoxy layers (using UMAT model).
Numerical results ( $η$ =0.001, 0^o loading angle)			Experimental results (De Vries, 2001)
C3D8R	C3D8	SC8R	Experimental results (De Vries, 2001)
463.7	467.1	458.7	446

Table 8. Net blunt notch strength (MPa) for different values of the viscosity parameter (using UMAT model).
Numerical results (C3D8R, 0^o loading angle)				Experimental results(De Vries, 2001)
$η$ =0.001	$η$ =0.0004	$η$ =0.00016	$η$ =0.000064	Experimental results(De Vries, 2001)
463.7	453.8	449.2	448.2	446

Table 9. Net blunt notch strength (MPa) for different values of the viscosity parameter in the fiber material (using multiscale model, viscosity in the matrix material $η_{m}$ =0.005).
Numerical results (C3D8R, 0^o loading angle)				Experimental results(De Vries, 2001)
$η$ =0.002	$η$ =0.001	$η$ =0.0005	$η$ =0.00025	Experimental results(De Vries, 2001)
442	432	423	415	446

Figures