- Email: [email protected]
FEM
Engineering Fracture Mechanics 72 (2005) 2549–2563 www.elsevier.com/locate/engfracmech
Numerical simulation study of fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch using combined BEM/FEM Hideki Sekine *, Bo Yan 1, Takeshi Yasuho Department of Aeronautics and Space Engineering, Tohoku University, 6-6-01 Aoba-yama, Aoba-ku, Sendai 980-8579, Japan Received 11 August 2003; received in revised form 28 August 2004; accepted 4 February 2005 Available online 13 June 2005
Abstract A combined boundary element method and finite element method (BEM/FEM) is employed to investigate the fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bonded fiber-reinforced polymer (FRP) composite patch. Numerical simulation of crack growth process of a cracked aluminum panel repaired with a FRP composite patch under uniaxial cyclic loading has been carried out. The curve of crack length on unpatched side of the cracked panel versus the number of cyclic loading is determined by the numerical simulation, and it agrees well with experimental data. Furthermore, the crack front profiles of the cracked panel during fatigue crack growth and the distributions of stress intensity factors along crack fronts are also numerically simulated. 2005 Elsevier Ltd. All rights reserved. Keywords: Cracked aluminum panel; FRP composite patch repair; Fatigue crack growth; Combined BEM/FEM
1. Introduction Repairs of cracked components in aerospace structures are becoming more and more important due to the requirement of operation safety. The repair methods based on adhesively bonded fiber-reinforced
*
1
Corresponding author. Tel.: +81 22 217 6982; fax: +81 22 217 6983. E-mail address: [email protected] (H. Sekine). On leave from Department of Engineering Mechanics, Chongqing University, Chongqing 400044, China.
0013-7944/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2005.02.007
2550
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
polymer (FRP) composite patches have been demonstrated to be very promising to these cracked structures [1,2]. FRP composite patches have the advantages of high ratios of stiffness and strength to weight, and are more structurally efficient and much less damaging to the repaired structures than fastened metallic patches. Although double-sided repair with FRP composite patches is more effective in reinforcement, single-sided repair plays a more important role because, in the most of practical repairs, it is difficult or even impossible to access both sides of the cracked structures that needed to be repaired. Fatigue crack growth behavior of cracked panels after being repaired decides the extension of fatigue life or service life of the repaired structures. Therefore, the evaluation of fatigue crack growth behavior of cracked panels repaired with a FRP composite patch becomes a focus in this research area. Some analytical and experimental studies have investigated the fatigue crack growth behavior of cracked metallic panels repaired with a single-sided FRP composite patch [2–5]. So far, the numerical studies for this problem have mainly been focused on the determination of stress intensity factors in the cracked panels. The first numerical model proposed by Ratwani [6] used two-dimensional finite elements to represent cracked panel and shear spring elements to represent adhesive layer by neglecting the influence of out-ofplane bending, which may lead to large errors in the most of practical situations. Several authors have used Mindlin plate finite elements to represent cracked panels in the numerical analyses [7–13]. With Mindlin plate finite element method (FEM), it is impossible to determine accurately the crack front profiles and the distributions of stress intensity factors along crack fronts. To evaluate accurately the fatigue crack growth behavior of the repaired cracked panels, it is necessary to determine the crack front profiles during crack growth and the distributions of stress intensity factors along crack fronts. Three-dimensional modeling of the cracked panels is a prerequisite for this purpose. The three-dimensional FEM has been used to analyze the single-sided repair problems by some authors [5,7,14]. However, these studies only dealt with the determination of stress intensity factors of cracked panels with a fixed crack length, no numerical simulation of the full process of fatigue crack growth has been performed. A combined boundary element method and finite element method (BEM/FEM) has been developed by Young [15]. With this method, the cracked panel was represented by three-dimensional boundary elements including traction singular quarter-point elements, and the FRP composite patch was represented by finite elements. The nodes on the attachment surfaces of the two portions, respectively represented by the boundary elements and the finite elements, are linked by means of springs which are used to represent the adhesive layer. Using this method, the stress intensity factors along crack front can be determined accurately and directly by boundary elements. Effective numerical investigation of fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch needs to simulate the fatigue crack growth process. However, except a few papers, e.g. [12], most of the works aforementioned focused on the determination of stress intensity factors of the repaired cracked panels at a certain crack length, but the crack growth process. In Ref. [12], stress intensity factors at the midplane of the cracked panel, which are calculated by Mindlin plate elements, have been used in the prediction of fatigue crack growth at the unpatched side, which may lead to some errors. This study, therefore, investigates numerically the fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch. Based on the facts that the stress intensity factors determined with the boundary element method (BEM) are more accurate than those with the FEM, when the number of the elements used in both models is the same, and the FEM is more powerful to analyze FRP composite laminates, the combined BEM/FEM with an effective method of automatic remeshing is adopted in the numerical simulation. Fatigue tests were performed to determine the material constants of fatigue crack growth and to provide the data for comparison with the results of numerical simulation. It has been found that the present numerical technique can appropriately evaluate the fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bonded FRP composite patch.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2551
2. Numerical simulation method using combined BEM/FEM The combined BEM/FEM proposed by Young [15], which is applied to analyze the stress intensity factors of cracked aluminum panels repaired with an adhesively bonded FRP composite patch, is briefly discussed, and an effective automatic remeshing method in the numerical simulation of the crack growth is then presented. A cracked aluminum panel repaired with a FRP composite patch is schematically shown in Fig. 1. 2.1. Boundary element method for cracked panels For an homogenous isotropic elastic body, the boundary element equation can be expressed as [16] Hub ¼ Gtb
ð1Þ
where ub and tb are vectors of nodal displacements and tractions respectively, and the matrices H and G are composed of integrals, which can be found in Ref. [16]. For the areas not adjacent to a crack front, the standard quadratic isoparametric elements can be used to obtain precise results. In the vicinity of crack front, the well-known traction singular quarter-point elements can be used to reflect the variation of displacements and the singularity of the tractions. If a set of orthogonal vectors n, q and p at a point on the crack front is defined as shown in Fig. 2, the three stress intensity factors KI, KII and KIII, corresponding to the opening, sliding and tearing modes, can be either expressed in terms of the modified tractions t* on the element ahead of the crack front as KI ¼
pffiffiffiffiffiffiffi 2plt n;
K II ¼
pffiffiffiffiffiffiffi 2plt q;
K III ¼
pffiffiffiffiffiffiffi 2plt p
Fig. 1. Cracked aluminum panel repaired with a FRP composite patch.
ð2Þ
2552
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563 n p
Crack Surfaces
q
Crack Front
Fig. 2. Direction vectors at crack front.
where l is the length of the side of the element ahead of the crack front, or displacement vectors as K I ¼ ku n;
K II ¼ ku q;
K III ¼ ð1 mÞku p
ð3Þ
8Dub ðl=4Þ Dub ðlÞ pffiffi 3 l
ð4Þ
in which 12 1 E 2pl ; k¼ 8 1 m2 lq
u ¼
where E and m are YoungÕs modulus and PoissonÕs ratio respectively, l is a vector from the nodal point on the crack front to the nodal point on the far edge of the element, and Dub(l/4) and Dub(l) are respectively the crack opening displacement vectors at the quarter-point distance and the far edge of the element from the crack front. It should be noted that at the intersection of crack front with a free edge, the tangential stress component will be finite, in which case the definitions of KI and KII must be multiplied by (1 m2) and KIII = 0. 2.2. Finite element method for FRP composite laminates The finite element formulation for FRP composite laminates is based on the extended Reissner variational principle as presented by Pagano [17]. Consider a laminate composed of NL layers lying in the (x, y) plane, with z denoting the third (normal) coordinate. The thickness of an individual layer with volume VM and surface AM is uniform. If the displacements and stresses are treated as independent variables, the equations of elasticity for the entire laminate can be deduced from the stationary point of the Hellinger– Reissner potential energy functional PHR as Z N L Z X 1 1 ~ ðui;j þ uj;i Þ rkl S ijkl rij dV M PHR ¼ ti ui dAM ð5Þ 2 2 VM AM M¼1 where ui is the displacement components, ui,j is the derivatives of displacements with respect to spatial coordinates, rkl the stress components, Sijkl the elastic compliance tensor and ~ti the tractions specified on surface AM. After introducing some generalized variables and discretizing the (x, y) domain of the laminate into quadratic isoparametric surface elements similar to the boundary elements, the potential energy functional for an individual element can be written as
T 1 T T T T T T PeHR ¼ Be Qe1 Ae þ Ce Qe2 Ae Be Qe4 Be þ 2Ce Qe5 Be þ Ce Qe6 Ce Ae f e þ Ce Qe3 De ð6Þ 2 where superscript e indicates elementary variables, Be is a nodal vector of stress averages and moments in an element and Ae is a nodal vector of displacement moments. The nodal vectors Ce and De are respectively the interlaminar stresses and displacements on the bottom and top surfaces, and Qek (k = 1–6) are coefficient matrices. Vector fe is equivalent nodal tractions on the edge surfaces of the layer corresponding to bound-
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2553
ary conditions. For the detailed expressions of the above vectors and matrices, one can refer to the work of Young [15]. Using Reissner variational principle for the entire domain, the global equation can be obtained in the form
QAA QTAC
QAC QCC
A C
¼
F QDC D
ð7Þ
where A and C are respectively the global nodal vectors of Ae and Ce in Eq. (6), and the global nodal vector B corresponding to Be is eliminated by means of element balance equations. The coefficient submatrices QAA, QAC, QCC and QDC are assemblies of the corresponding element submatrices, which are dependent on Qek (i = 1–6), and F and D are the global vectors of the boundary tractions fe on edge surfaces and displacements De on surfaces of a layer respectively. It should be noted that the interfacial displacement components of D are either absent from the right side of the governing equation (7) or are specified at some interfaces. This means that the components of D only on the external surfaces of the composite laminate are kept in Eq. (7), which are unknowns or specified by boundary conditions, in the case of no interfacial values being specified. 2.3. Linking of boundary elements and finite elements by means of springs In practical situations, adhesive layer is usually much thinner and more flexible than the cracked aluminum panels and FRP composite patches. It is reasonable to simplify the adhesive layer as springs possessing shear stresses r13 and r23 and normal stress r33 only. For the modeling of adhesive layer with linear springs, the relationship between the tractions Ts in springs and the relative displacements Us of the both ends of the springs is Us ¼ qTs
ð8Þ n
where q is a block diagonal matrix assembled from the nodal 3 · 3 adhesive flexibility matrices q , which is diagonal and depends on the thickness and elastic constants of the adhesive layer. The nodal parameters C and D in Eq. (7) of the finite element portion of the model can be subdivided into those on the attachment surface denoted by Tc and Uc respectively, and the rest corresponding to specified boundary conditions or interlayer tractions, C 0 and D 0 . In the rest of this paper, the superscript c indicates the variables of the nodes on attachment surface. It is noted that the submatrix QDC in Eq. (7) is then subdivided into two submatrices QDT and QUT corresponding to D 0 and Uc respectively. Grouping together all the unknowns A and C 0 into the set A 0 for convenience, and collecting the known values of A and C 0 along with the loading terms F and QDTD 0 into the vectors FA and FT, the equation system (7) may be rearranged into the form QA0 A0 A0 þ QA0 T Tc ¼ FA QTA0 T A0 þ QTT Tc ¼ FT QUT Uc
ð9Þ
In order to keep only the variables on the attachment surface, the vector A 0 in Eq. (9) can be eliminated, and the equations are then simplified to Uc ¼ Qc Tc þ F c
ð10Þ
where Q and F are respectively a matrix and a vector depending on QUT, QA 0 A 0 , QA 0 T, QTT, FA and FT. It is noted that the matrix Qc is generally asymmetric and fully populated. For the boundary element portion of the model, the boundary element system of Eq. (1) can be rearranged into
2554
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
H c uc G c t c þ H o uo G o t o ¼ R
ð11Þ
where the superscript o indicates the variables on the remainder surface of the boundary element portion other than the attachment surface, and R includes all the terms relating to the specified boundary conditions. This modified boundary element system contains more unknowns than those in regular boundary element equations due to the presence of both uc and tc as unknowns on the attachment surface. Now the nodes on the attachment surfaces of the boundary element portion and the finite element portion are linked by means of the springs representing adhesive layer, and the conditions of continuity and compatibility can be expressed as tc ¼ Tc c
c
u ¼U þU
ð12Þ s
ð13Þ s
c
Substituting Eq. (8) into Eq. (13) and considering T = T , Eq. (13) becomes uc ¼ Uc þ qTc
ð14Þ
Combination of Eqs. (10)–(12) and (14) leads to the following system equation: Ho uo Go to ½Gc þ Hc ðQc þ qÞ tc ¼ R Hc F
ð15Þ
which is applied to analyze the stress intensity factors of the cracked panels repaired with a FRP composite patch. 2.4. Automatic remeshing An effective method is presented to automatically remesh the model in numerical simulation of the fatigue crack growth in cracked panels repaired with a FRP composite patch. In the numerical simulation, both the boundary element mesh of the cracked panel and the finite element mesh of the composite patch are remeshed in each increment step of the number of cyclic loading. A typical boundary element mesh of the cross section of the cracked panel, on which the crack surface is on the left of the crack front, is shown in Fig. 3. On this cross section, only the rectangular zone a 0 bcd 0 is remeshed during the numerical simulation of the crack growth. In each increment step of the number of cyclic loading, the increment of crack length at each node on the crack front is determined based on Paris law, which is discussed in detail in Section 3.1. Using the increment of crack length and the normal direction of the current crack front at each node, the new position of each node on the crack front can be determined. The new crack front profile obtained by smoothing the curve determined with the new coordinates of the nodes on the crack front is then discretized by the nodes uniformly located along the thickness direction of the cracked panel. To ensure the accuracy of the numerical simulation, the sizes of the two column elements
Traction Singular Quarter-Point Element b
a
b’
a’
c
d
c’
d’
Crack Front
Fig. 3. Boundary element mesh of cross section with crack surface.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2555
on each side of the crack front, i.e., the elements in the zone adc 0 b 0 , are smaller than those of the others and keep the same during the simulation of the fatigue crack growth. It is noted that the coordinates of the nodes along the lines bc and a 0 d 0 , as shown in Fig. 3, keep the same during the numerical simulation of the crack growth, and the positions of the nodes along edges ad and b 0 c 0 change with the growth of the fatigue crack. In addition, the new positions of the other nodes in zones abcd and a 0 b 0 c 0 d 0 are automatically generated at each increment step of the number of cyclic loading by means of isoparametric transformation, and a bias factor is used to ensure the elements close to the crack front being smaller than those far away from the crack front. On the top and bottom surfaces of the cracked panel, the boundary elements which include the nodes along the edges ba 0 and cd 0 are modified in each step. In those elements, the new positions of the midpoint nodes of each element are located at the midpoints of the four edges connecting the four new corner nodes of the element. It is noted that the small elements around the crack tips on the top and bottom surfaces of the cracked panel keep their shapes and move with the crack tips in the numerical simulation. Meanwhile, the boundary element mesh of the other part of the cracked panel keeps the same in each step. After remeshing the boundary element portion of the model in an increment step of the number of cyclic loading, the finite element mesh of the composite patch is updated simply by replacing the x and y coordinates of the finite element nodes with the new x and y coordinates of the corresponding boundary element nodes on the attachment surface of the cracked panel. Some numerical examples have been used to demonstrate that this remeshing method is effective for this typical problem. It is noted that if the shapes of the boundary elements on the top and bottom surfaces of the cracked panel, which include the nodes along the edges ba 0 and cd 0 as shown in Fig. 3, and the corresponding finite elements of the composite patch become irregular due to a long distant propagation of the crack, a manual remeshing is needed at this step.
3. Numerical models and experiments 3.1. Numerical models The fatigue crack growth behavior of a cracked aluminum panel repaired with a circular FRP composite patch, as shown in Fig. 1, is numerically simulated. The material of the cracked panel is 2024-T3 aluminum and the adhesive is AF-163-2K. The composite patch is a ½0 3 = 45 =0 2 s glass/epoxy laminate, where the fiber direction of 0 is parallel to the longitudinal direction of the panel. The dimensions of the cracked aluminum panel, FRP composite patch and adhesive layer are listed in Table 1, and the corresponding mechanical properties are shown in Table 2. The cracked aluminum panel is subject to a cyclic stress with amplitude Dr0 = 59.6 MPa, and the applied maximum stress is rmax = 79.4 MPa and the stress ratio R = 0.25. The initial crack length on the surface of this cracked panel repaired with composite patch is 2a = 21.30 mm. In the numerical simulation, the taper of the patch, as shown in Fig. 1, is ignored and its diameter is assumed to be D. The adhesive layer is simplified as linear springs as discussed in Section 2.3. To reduce
Table 1 Dimensions of cracked panel with a FRP composite patch Components
Dimensions
Cracked aluminum panel FRP composite patch Adhesive layer
L = 248 mm, W = 118 mm, tp = 2 mm D = 50 mm, tl = 3.2 mm ta = 0.1 mm
2556
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
Table 2 Mechanical properties of cracked panel with a FRP composite patch Materials
Mechanical properties
2024-T3 aluminum Lamina of composite AF-163-2K adhesive
E = 72.39 GPa, m = 0.33 E11 = 44.12 GPa, E22 = 9.65 GPa, G12 = 4.13 GPa, m12 = 0.28 Ea = 1.18 GPa, ma = 0.34
the computational cost, an equivalent transversely isotropic single layer is employed to represent the original FRP composite patch. In order to obtain the effective material properties of this equivalent single layer, an optimization problem has been set up, in which the objective function is chosen as the sum of the leastsquares of differences between the in-plane stiffness and bending stiffness of this single layer and those of the original composite laminate, and the design variables are the material constants. Through the minimization of this objective function, the material constants of the FRP composite patch represented by the equivalent single layer are obtained as E11 ¼ 37.35 GPa;
E22 ¼ 11.38 GPa;
G12 ¼ 5.97 GPa;
m12 ¼ 0.38
where subscripts 1 and 2 are respectively parallel to the longitudinal and transverse directions of the panel. Owing to the symmetry of the problem, only one-fourth of the model is discretized. The boundary element mesh and finite element mesh are shown in Fig. 4. For the boundary element portion of the model, at the symmetric surface x = 0, the nodal displacements ux = 0, and at the symmetric surface y = 0, except for the crack surface, uy = 0. In the finite element portion, the symmetric boundary conditions are set at the two symmetric edges of the composite patch and the displacement in z direction at the top of the node at the center of the circular patch is set to zero. The average stress is applied on the loading end of the cracked panel. It is noted that the dimension in thickness direction of the schematic model in Fig. 4 is amplified by three times for clarity of view. It is important to investigate how many boundary element layers in the thickness direction of the cracked panel shown in Fig. 4, around the crack front area, are adequate to get accurate stress intensity factors along the crack front. Several mesh models of the unpached cracked panel with different numbers of boundy
FRP Composite Patch
z
Aluminum Panel
Crack Surface Crack Front
x
Fig. 4. BEM/FEM mesh of cracked aluminum panel repaired with a FRP composite patch.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2557
Table 3 Stress intensity factors of cracked aluminum panel with a FRP composite patch z (mm)
KI (MPa m1/2) with BEM/FEM
KI (MPa m1/2) with FEM
Error (%)
0.40 0.80 1.20 1.60
12.26 10.98 9.46 7.72
12.60 11.09 9.52 7.84
1.72 1.00 0.63 1.59
ary element layers in the thickness direction of the cracked panel are analyzed by means of the BEM. After comparison of these results, it is concluded that five element layers in the thickness direction of the cracked panel are sufficient to ensure the accuracy of the numerical simulation in this case. To demonstrate the validity of the combined BEM/FEM and the mesh of the model, as shown in Fig. 4, a model with a straight crack of the length a = 15 mm is analyzed by means of the combined BEM/FEM and the FEM. The mesh of the combined BEM/FEM model is similar to that as in Fig. 4. The mesh of the FEM model is much fine, which is composed of 61,046 nodes and 13,292 three-dimensional 20-node isoparametric elements, and the singular quarter-point elements are used to model the crack front area. The stress intensity factors at the points corresponding to the boundary element corner nodes along the crack front, except for the cross points of the crack front with the top and bottom surfaces of the cracked panel, obtained with these two methods are compared in Table 3. It is shown that the difference between the stress intensity factors obtained with the combined BEM/FEM and the FEM are small. The agreement of the stress intensity factors at the points on the crack front demonstrates the validity of the combined BEM/ FEM and the mesh model used in the numerical simulation of fatigue crack growth in the cracked aluminum panel repaired with a FRP composite patch. In the numerical simulation of the fatigue crack growth, the stress intensity factor KI max of the cracked panel with a certain crack front profile under the action of maximum stress rmax is firstly determined by means of the combined BEM/FEM. Because the stress intensity factors at the cross points of the crack front with the top and bottom surfaces of the cracked panel cannot be obtained accurately, the stress intensity factors at these points are modified by smoothing with the values at other nodes along the crack front. Taking account of Eqs. (3) and (4), the stress intensity factor range DKI corresponding to the stress range Dr0 can be determined by DK I ¼ ð1 RÞK I max
ð16Þ
The increment of crack length at each node along crack front, corresponding to an increment of the number of cyclic loading, can then be determined based on Paris law da ¼ CðDK I Þm dN
ð17Þ
where da/dN is the crack growth rate, and C and m are the material constants of the aluminum panels. In addition, because the influence of residual stress generated in the process of patching FRP composite laminates on stress intensity factor range is not great [1], the residual stress is neglected in the numerical simulation. The disbond effect of patch during fatigue crack growth process is also ignored in the numerical simulation for simplification. 3.2. Experiments To determine the material constants of fatigue crack growth of 2024-T3 aluminum, fatigue test of an unpatched center-cracked aluminum panel was performed. In this test, the cracked aluminum panel was subject to a sinusoidal stress with amplitude Dr0, a stress ratio of 0.25 and a cyclic frequency of 0.5 Hz
2558
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
Table 4 Dimensions and loading of unpatched cracked aluminum panel specimen L (mm)
W (mm)
tp (mm)
2a (mm)
Dr0 (MPa)
248
118
5.0
60.0
23.8
40
a (mm)
35
30
Numerical result Experimental data
25
20 0
1
2
3
4
[×104]
N (cycles) Fig. 5. Crack length of unpatched cracked aluminum panel versus number of cyclic loading.
in the environmental temperature of 27 C. The dimensions of the specimen, the initial crack length 2a and the applied cyclic stress Dr0 are listed in Table 4. The experimental data of the crack length versus the number of cyclic loading are shown in Fig. 5. Fatigue test of the cracked aluminum panel repaired with a circular FRP composite patch, as shown in Fig. 1, was then performed. The materials and dimensions of the cracked panel, adhesive and composite patch, as well as the applied stress are the same as those of the model discussed in Section 3.1. The other test conditions are the same as those in the test of the unpatched cracked aluminum panel mentioned above. Before the cracked aluminum panel is repaired, it was applied some cycles of loading until the crack length on the surface arrives at 2a = 21.30 mm. The cracked aluminum panel was then repaired with a circular FRP composite patch and was continuously applied the same cyclic loading until fatigue failure. The experimental data of the crack length on the unpatched side of the cracked panel repaired with a FRP composite patch versus the number of cyclic loading are shown in Fig. 6. For the purpose to investigate the fatigue crack growth behavior, the specimens of the unpatched and repaired aluminum cracked panels were fractured after some cycles of loading, and the crack front profiles of the specimens were then observed. The micrographs of typical examples of the crack front profiles corresponding to these two cases are respectively shown in Figs. 7 and 8.
4. Numerical results and discussion Fitting the results of the numerical simulations of fatigue crack growth of the unpatched cracked panel with the BEM to the experiment data gives the material constants of 2024-T3 aluminum C = 7.04 · 107
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2559
25
∆σ0=80.0 MPa 20
∆σ0=59.6 MPa
a (mm)
15
∆σ0=40.0 MPa
10
Present numerical model Simplified numerical model Experimental data
5
0 0
1
2
3
4
5
6
7
8
9
10 [×104]
N (cycles) Fig. 6. Crack length on unpatched side of cracked aluminum panel versus number of cyclic loading.
Fig. 7. Crack front of unpatched specimen (material: 2024-T3 aluminum; thickness: 5 mm; number of cyclic loading: 46448).
(MPa m1/2)2.48 mm/cycle and m = 2.48. In the numerical simulations, the mesh of the boundary element model is similar with those of the boundary element portion of the combined BEM/FEM model as shown in Fig. 4. Besides the boundary conditions as set in the boundary element portion of the model in Fig. 4, the displacement in z direction of the node located at the origin of the coordinate system is set to zero in this BEM model. The variation of crack length a versus the number of cyclic loading N determined by the numerical simulation is illustrated in Fig. 5, and it agrees well with the experimental data. In addition, the crack front profiles and the distributions of stress intensity factors of the unpatched specimen under the maximum stress rmax = 31.76 MPa at some cycles of loading during the fatigue crack growth simulated with the BEM are shown in Fig. 9. It is obvious that with the change of the distribution of stress intensity factors along crack front, the crack front is gradually curved from the initial straight line and propagates nearly parallel after some cycles of loading. The shapes of the crack fronts in the process of stable propagation determined by numerical simulation as in Fig. 9(a) is similar to that observed in experiment as in Fig. 7.
2560
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
z (mm)
Fig. 8. Crack front of specimen repaired with a FRP composite patch (material: 2024-T3 aluminum; thickness: 2 mm; number of cyclic loading: 61,600).
5 4 3 2 1 0 22
A
24
26
E
D
C
B
28
30
32
F
34
36
G
38
40
42
x (mm)
(a) 20.0
G F
15.0
KI (MPa m1/2)
E D 10.0
A
5.0
B
C
0.0 0.0
(b)
1.0
2.0
3.0
4.0
5.0
z (mm)
Fig. 9. Crack fronts and stress intensity factors in the case of unpatched cracked aluminum panel (number of cyclic loading: A: 0; B: 1000; C: 8000; D: 25,000; E: 49,000; F: 66,000; G: 78,000). (a) Crack fronts and (b) stress intensity factors.
Using the material constants, the combined BEM/FEM is employed to numerically simulate the fatigue crack growth behavior of the cracked aluminum panel repaired with a FRP composite patch. It is noted that the fatigue crack growth process of the cracked aluminum panel before being repaired is analyzed with
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2561
the BEM, and the final crack front profile in this BEM simulation is adopted as the initial profile of the repaired cracked panel. The variation of crack length on the unpatched side of the cracked aluminum panel repaired with a FRP composite patch, in the case of Dr0 = 59.6 MPa, versus the number of cyclic loading, which is determined by the numerical simulation, is demonstrated in Fig. 6. It is shown that the curve determined with the numerical method agrees well with the experiment data. The numerically simulated crack front profiles and distributions of stress intensity factors of the cracked aluminum panel under the peak stress rmax = 79.4 MPa are illustrated in Fig. 10. After the cracked aluminum panel is repaired, the deformation of the cracked panel is constrained by the FRP composite patch, so that the stress intensity factors along crack front decrease. The constraint of the FRP composite patch is implemented through the shear force transferred by adhesive layer. At the very beginning of fatigue crack growth after its reinitiating, the cross point of crack front with the unpatched surface of the cracked panel grows quickly due to the larger stress intensity factor, comparing with other points along the crack front. Meanwhile, the stress intensity factor at the cross point decreases and those at the other points increase as the number of cyclic loading increases. Therefore, the growth rate of crack on the unpatched side decreases slightly during the initial cyclic loading. By comparing Figs. 8 and 10(a), it is found that the crack front profiles determined with the combined BEM/FEM numerical simulation are similar to that observed in experiment except the portion close to the unpatched side of the cracked panel. As aforementioned, most of the previous works such as Ref. [12] did not consider the variation of stress intensity factors along curved crack front. However, as illustrated in Figs. 8 and 10, the crack front profiles of
z (mm)
2.0 1.5
B
1.0
A
0.5 0.0 10
D
C 11
12
13
14
F
E 15
16
17
G
18
19
20
x (mm)
(a) 15.0
F
10.0
E
D
1/2
KI (MPa m )
G
C B
5.0
A
0.0 0.0
(b)
0.5
1.0
1.5
2.0
z (mm)
Fig. 10. Crack fronts and stress intensity factors of cracked aluminum panel repaired with a FRP composite patch (number of cyclic loading: A: 0; B: 4000; C: 8000; D: 30,000; E: 45,000; F: 66,000; G: 82,000). (a) Crack fronts and (b) stress intensity factors.
2562
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
cracked panel repaired with a FRP composite patch are not straight lines and the distributions of the stress intensity factors along the crack front are nonlinear during the fatigue crack growth. Therefore, it is important to consider these two factors in the numerical simulation of the fatigue crack growth behavior of the cracked panel repaired with a FRP composite patch. Using the average of stress intensity factors along the crack fronts, except the values at cross points of the crack front with the top and bottom surfaces of the cracked panel, and assuming that the crack front keep straight and perpendicular to the top and bottom surfaces of the cracked panel and propagate parallel during the crack growth, the fatigue crack growth of the same model as shown in Fig. 4, under the cyclic stress Dr0 = 59.6 MPa, is numerically simulated by means of the combined BEM/FEM. For the convenience of description, this model is referred as the simplified numerical model. The curve of the crack length versus the number of cyclic loading of the simplified numerical model is shown in Fig. 6 with a dash line. It is obvious that the curve of the simplified numerical model is close to the experimental data at the initial stage, but the error becomes larger and larger with the increase of the number of cyclic loading. Although the results of the simplified numerical model are closer to the experimental data, in the initial stage, than those of the present numerical model in which the curved crack front and the nonlinear distribution of the stress intensity factors along the crack front are taken into account, the results of the later are closer to the experimental data in the whole stage of the fatigue crack growth. To investigate the fatigue crack growth behavior of the cracked panel repaired with a FRP composite patch under two different cyclic stresses with amplitudes Dr0 = 40.0 MPa and 80.0 MPa, the numerical simulations are carried out with the combined BEM/FEM. The variations of the crack length on the unpatched side of the cracked aluminum panel versus the number of cyclic loading in these two cases are shown in Fig. 6. It is illustrated that as the increase of the amplitude of the cyclic loading, the crack front propagates much faster.
5. Conclusions Fatigue crack growth behavior of a cracked aluminum panel repaired with a FRP composite patch under uniaxial cyclic loading is numerically investigated by the combined BEM/FEM. The curve of fatigue crack length on the unpatched side of the cracked panel versus the number of cyclic loading, which is determined by the numerical simulation, agrees well with experimental data. The crack front profiles and the distributions of stress intensity factors along the crack fronts during fatigue crack growth are numerically simulated. All these results demonstrate that the numerical method used in this paper is efficient in the evaluation of fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bounded FRP composite patch.
Acknowledgements This work was partly supported by the Grand-in-Aid for Scientific Research (A) (2) No. 14205138 and (B) (2) No. 13555264 from the Ministry of Education, Culture, Sports, Science and Technology of Japan to HS. The authors would like to thank Mr. Takamitsu Sato for his assistance to conducting the experiments.
References [1] Baker AA. Repair efficiency in fatigue-cracked aluminum components. Fatigue Fract Engng Mater Struct 1993;16(7):753–65. [2] Baker AA, Rose LRF, Jones R. Advances in the bonded composite repairs of metallic aircraft structures, vols. 1&2. London: Elsevier Science Publications; 2002.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2563
[3] Denney JJ, Mall S. Characteristics of disbanded effects on fatigue crack growth behavior in aluminum plate with bonded composite patch. Engng Fract Mech 1997;57(5):507–25. [4] Schubbe JJ, Mall S. Investigation of a cracked thick aluminum panel repaired with a bonded composite patch. Engng Fract Mech 1999;63(3):305–23. [5] Seo DC, Lee JJ. Fatigue crack growth behavior of cracked aluminum plate repaired with composite patch. Compos Struct 2002;57(1–4):323–30. [6] Ratwani MM. Analysis of cracked adhesively bonded laminated structures. AIAA J 1979;17(9):988–94. [7] Sun CT, Klug J, Arendt C. Analysis of cracked aluminum plates repaired with bonded composite panels. AIAA J 1996;34(2):369–74. [8] Naboulsi S, Mall S. Modeling of a cracked metallic structure with bonded composite patch using the three-layer technique. Compos Struct 1996;35(3):295–308. [9] Naboulis S, Mall S. Nonlinear analysis of bonded composite patch repair of cracked aluminum panels. Compos Struct 1998; 41(3–4):303–13. [10] Schubbe JJ, Mall S. Modeling of cracked thick metallic structure with bonded composite patch repair using three-layer technique. Compos Struct 1999;45(3):185–93. [11] Naboulsi S, Mall S. Fatigue crack growth analysis of adhesively repaired panel using perfectly and imperfectly composite patches. Theor Appl Fract Mech 1997;28(1):13–28. [12] Shibuya Y, Fujimoto S, Aoki D, Sato M, Shirahata H, Fukunaga H, et al. Evaluation of crack growth in cracked aluminum panels repaired with a bonded composite patch under cyclic loading. Adv Compos Mater 2001;10(4):287–97. [13] Kam TY, Chu KH, Tsai YC. Fatigue of cracked plates repaired with single-sided composite patch. AIAA J 1998;36(4):645–50. [14] Chue CH, Chang LC, Tsai JS. Bonded repair of a plate with inclined central crack under biaxial loading. Compos Struct 1994;28(1):39–45. [15] Young A. Three-dimensional analysis of a patched cracked solid using a combined boundary element and finite element model. DRP 26C-R02-Milestone 2 (Item 410), Working Paper MS4-93-WP-24, Materials and Structures Department, Defense Research Agency, Hampshire, 1993. [16] Brebbia CA. The boundary element method for engineering. London: Pentech Press; 1978. [17] Pagano NJ. Stress field in composite laminates. Int J Solids Struct 1978;14(5):385–400.
Numerical simulation study of fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch using combined BEM/FEM Hideki Sekine *, Bo Yan 1, Takeshi Yasuho Department of Aeronautics and Space Engineering, Tohoku University, 6-6-01 Aoba-yama, Aoba-ku, Sendai 980-8579, Japan Received 11 August 2003; received in revised form 28 August 2004; accepted 4 February 2005 Available online 13 June 2005
Abstract A combined boundary element method and finite element method (BEM/FEM) is employed to investigate the fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bonded fiber-reinforced polymer (FRP) composite patch. Numerical simulation of crack growth process of a cracked aluminum panel repaired with a FRP composite patch under uniaxial cyclic loading has been carried out. The curve of crack length on unpatched side of the cracked panel versus the number of cyclic loading is determined by the numerical simulation, and it agrees well with experimental data. Furthermore, the crack front profiles of the cracked panel during fatigue crack growth and the distributions of stress intensity factors along crack fronts are also numerically simulated. 2005 Elsevier Ltd. All rights reserved. Keywords: Cracked aluminum panel; FRP composite patch repair; Fatigue crack growth; Combined BEM/FEM
1. Introduction Repairs of cracked components in aerospace structures are becoming more and more important due to the requirement of operation safety. The repair methods based on adhesively bonded fiber-reinforced
*
1
Corresponding author. Tel.: +81 22 217 6982; fax: +81 22 217 6983. E-mail address: [email protected] (H. Sekine). On leave from Department of Engineering Mechanics, Chongqing University, Chongqing 400044, China.
0013-7944/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2005.02.007
2550
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
polymer (FRP) composite patches have been demonstrated to be very promising to these cracked structures [1,2]. FRP composite patches have the advantages of high ratios of stiffness and strength to weight, and are more structurally efficient and much less damaging to the repaired structures than fastened metallic patches. Although double-sided repair with FRP composite patches is more effective in reinforcement, single-sided repair plays a more important role because, in the most of practical repairs, it is difficult or even impossible to access both sides of the cracked structures that needed to be repaired. Fatigue crack growth behavior of cracked panels after being repaired decides the extension of fatigue life or service life of the repaired structures. Therefore, the evaluation of fatigue crack growth behavior of cracked panels repaired with a FRP composite patch becomes a focus in this research area. Some analytical and experimental studies have investigated the fatigue crack growth behavior of cracked metallic panels repaired with a single-sided FRP composite patch [2–5]. So far, the numerical studies for this problem have mainly been focused on the determination of stress intensity factors in the cracked panels. The first numerical model proposed by Ratwani [6] used two-dimensional finite elements to represent cracked panel and shear spring elements to represent adhesive layer by neglecting the influence of out-ofplane bending, which may lead to large errors in the most of practical situations. Several authors have used Mindlin plate finite elements to represent cracked panels in the numerical analyses [7–13]. With Mindlin plate finite element method (FEM), it is impossible to determine accurately the crack front profiles and the distributions of stress intensity factors along crack fronts. To evaluate accurately the fatigue crack growth behavior of the repaired cracked panels, it is necessary to determine the crack front profiles during crack growth and the distributions of stress intensity factors along crack fronts. Three-dimensional modeling of the cracked panels is a prerequisite for this purpose. The three-dimensional FEM has been used to analyze the single-sided repair problems by some authors [5,7,14]. However, these studies only dealt with the determination of stress intensity factors of cracked panels with a fixed crack length, no numerical simulation of the full process of fatigue crack growth has been performed. A combined boundary element method and finite element method (BEM/FEM) has been developed by Young [15]. With this method, the cracked panel was represented by three-dimensional boundary elements including traction singular quarter-point elements, and the FRP composite patch was represented by finite elements. The nodes on the attachment surfaces of the two portions, respectively represented by the boundary elements and the finite elements, are linked by means of springs which are used to represent the adhesive layer. Using this method, the stress intensity factors along crack front can be determined accurately and directly by boundary elements. Effective numerical investigation of fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch needs to simulate the fatigue crack growth process. However, except a few papers, e.g. [12], most of the works aforementioned focused on the determination of stress intensity factors of the repaired cracked panels at a certain crack length, but the crack growth process. In Ref. [12], stress intensity factors at the midplane of the cracked panel, which are calculated by Mindlin plate elements, have been used in the prediction of fatigue crack growth at the unpatched side, which may lead to some errors. This study, therefore, investigates numerically the fatigue crack growth behavior of cracked aluminum panels repaired with a FRP composite patch. Based on the facts that the stress intensity factors determined with the boundary element method (BEM) are more accurate than those with the FEM, when the number of the elements used in both models is the same, and the FEM is more powerful to analyze FRP composite laminates, the combined BEM/FEM with an effective method of automatic remeshing is adopted in the numerical simulation. Fatigue tests were performed to determine the material constants of fatigue crack growth and to provide the data for comparison with the results of numerical simulation. It has been found that the present numerical technique can appropriately evaluate the fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bonded FRP composite patch.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2551
2. Numerical simulation method using combined BEM/FEM The combined BEM/FEM proposed by Young [15], which is applied to analyze the stress intensity factors of cracked aluminum panels repaired with an adhesively bonded FRP composite patch, is briefly discussed, and an effective automatic remeshing method in the numerical simulation of the crack growth is then presented. A cracked aluminum panel repaired with a FRP composite patch is schematically shown in Fig. 1. 2.1. Boundary element method for cracked panels For an homogenous isotropic elastic body, the boundary element equation can be expressed as [16] Hub ¼ Gtb
ð1Þ
where ub and tb are vectors of nodal displacements and tractions respectively, and the matrices H and G are composed of integrals, which can be found in Ref. [16]. For the areas not adjacent to a crack front, the standard quadratic isoparametric elements can be used to obtain precise results. In the vicinity of crack front, the well-known traction singular quarter-point elements can be used to reflect the variation of displacements and the singularity of the tractions. If a set of orthogonal vectors n, q and p at a point on the crack front is defined as shown in Fig. 2, the three stress intensity factors KI, KII and KIII, corresponding to the opening, sliding and tearing modes, can be either expressed in terms of the modified tractions t* on the element ahead of the crack front as KI ¼
pffiffiffiffiffiffiffi 2plt n;
K II ¼
pffiffiffiffiffiffiffi 2plt q;
K III ¼
pffiffiffiffiffiffiffi 2plt p
Fig. 1. Cracked aluminum panel repaired with a FRP composite patch.
ð2Þ
2552
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563 n p
Crack Surfaces
q
Crack Front
Fig. 2. Direction vectors at crack front.
where l is the length of the side of the element ahead of the crack front, or displacement vectors as K I ¼ ku n;
K II ¼ ku q;
K III ¼ ð1 mÞku p
ð3Þ
8Dub ðl=4Þ Dub ðlÞ pffiffi 3 l
ð4Þ
in which 12 1 E 2pl ; k¼ 8 1 m2 lq
u ¼
where E and m are YoungÕs modulus and PoissonÕs ratio respectively, l is a vector from the nodal point on the crack front to the nodal point on the far edge of the element, and Dub(l/4) and Dub(l) are respectively the crack opening displacement vectors at the quarter-point distance and the far edge of the element from the crack front. It should be noted that at the intersection of crack front with a free edge, the tangential stress component will be finite, in which case the definitions of KI and KII must be multiplied by (1 m2) and KIII = 0. 2.2. Finite element method for FRP composite laminates The finite element formulation for FRP composite laminates is based on the extended Reissner variational principle as presented by Pagano [17]. Consider a laminate composed of NL layers lying in the (x, y) plane, with z denoting the third (normal) coordinate. The thickness of an individual layer with volume VM and surface AM is uniform. If the displacements and stresses are treated as independent variables, the equations of elasticity for the entire laminate can be deduced from the stationary point of the Hellinger– Reissner potential energy functional PHR as Z N L Z X 1 1 ~ ðui;j þ uj;i Þ rkl S ijkl rij dV M PHR ¼ ti ui dAM ð5Þ 2 2 VM AM M¼1 where ui is the displacement components, ui,j is the derivatives of displacements with respect to spatial coordinates, rkl the stress components, Sijkl the elastic compliance tensor and ~ti the tractions specified on surface AM. After introducing some generalized variables and discretizing the (x, y) domain of the laminate into quadratic isoparametric surface elements similar to the boundary elements, the potential energy functional for an individual element can be written as
T 1 T T T T T T PeHR ¼ Be Qe1 Ae þ Ce Qe2 Ae Be Qe4 Be þ 2Ce Qe5 Be þ Ce Qe6 Ce Ae f e þ Ce Qe3 De ð6Þ 2 where superscript e indicates elementary variables, Be is a nodal vector of stress averages and moments in an element and Ae is a nodal vector of displacement moments. The nodal vectors Ce and De are respectively the interlaminar stresses and displacements on the bottom and top surfaces, and Qek (k = 1–6) are coefficient matrices. Vector fe is equivalent nodal tractions on the edge surfaces of the layer corresponding to bound-
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2553
ary conditions. For the detailed expressions of the above vectors and matrices, one can refer to the work of Young [15]. Using Reissner variational principle for the entire domain, the global equation can be obtained in the form
QAA QTAC
QAC QCC
A C
¼
F QDC D
ð7Þ
where A and C are respectively the global nodal vectors of Ae and Ce in Eq. (6), and the global nodal vector B corresponding to Be is eliminated by means of element balance equations. The coefficient submatrices QAA, QAC, QCC and QDC are assemblies of the corresponding element submatrices, which are dependent on Qek (i = 1–6), and F and D are the global vectors of the boundary tractions fe on edge surfaces and displacements De on surfaces of a layer respectively. It should be noted that the interfacial displacement components of D are either absent from the right side of the governing equation (7) or are specified at some interfaces. This means that the components of D only on the external surfaces of the composite laminate are kept in Eq. (7), which are unknowns or specified by boundary conditions, in the case of no interfacial values being specified. 2.3. Linking of boundary elements and finite elements by means of springs In practical situations, adhesive layer is usually much thinner and more flexible than the cracked aluminum panels and FRP composite patches. It is reasonable to simplify the adhesive layer as springs possessing shear stresses r13 and r23 and normal stress r33 only. For the modeling of adhesive layer with linear springs, the relationship between the tractions Ts in springs and the relative displacements Us of the both ends of the springs is Us ¼ qTs
ð8Þ n
where q is a block diagonal matrix assembled from the nodal 3 · 3 adhesive flexibility matrices q , which is diagonal and depends on the thickness and elastic constants of the adhesive layer. The nodal parameters C and D in Eq. (7) of the finite element portion of the model can be subdivided into those on the attachment surface denoted by Tc and Uc respectively, and the rest corresponding to specified boundary conditions or interlayer tractions, C 0 and D 0 . In the rest of this paper, the superscript c indicates the variables of the nodes on attachment surface. It is noted that the submatrix QDC in Eq. (7) is then subdivided into two submatrices QDT and QUT corresponding to D 0 and Uc respectively. Grouping together all the unknowns A and C 0 into the set A 0 for convenience, and collecting the known values of A and C 0 along with the loading terms F and QDTD 0 into the vectors FA and FT, the equation system (7) may be rearranged into the form QA0 A0 A0 þ QA0 T Tc ¼ FA QTA0 T A0 þ QTT Tc ¼ FT QUT Uc
ð9Þ
In order to keep only the variables on the attachment surface, the vector A 0 in Eq. (9) can be eliminated, and the equations are then simplified to Uc ¼ Qc Tc þ F c
ð10Þ
where Q and F are respectively a matrix and a vector depending on QUT, QA 0 A 0 , QA 0 T, QTT, FA and FT. It is noted that the matrix Qc is generally asymmetric and fully populated. For the boundary element portion of the model, the boundary element system of Eq. (1) can be rearranged into
2554
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
H c uc G c t c þ H o uo G o t o ¼ R
ð11Þ
where the superscript o indicates the variables on the remainder surface of the boundary element portion other than the attachment surface, and R includes all the terms relating to the specified boundary conditions. This modified boundary element system contains more unknowns than those in regular boundary element equations due to the presence of both uc and tc as unknowns on the attachment surface. Now the nodes on the attachment surfaces of the boundary element portion and the finite element portion are linked by means of the springs representing adhesive layer, and the conditions of continuity and compatibility can be expressed as tc ¼ Tc c
c
u ¼U þU
ð12Þ s
ð13Þ s
c
Substituting Eq. (8) into Eq. (13) and considering T = T , Eq. (13) becomes uc ¼ Uc þ qTc
ð14Þ
Combination of Eqs. (10)–(12) and (14) leads to the following system equation: Ho uo Go to ½Gc þ Hc ðQc þ qÞ tc ¼ R Hc F
ð15Þ
which is applied to analyze the stress intensity factors of the cracked panels repaired with a FRP composite patch. 2.4. Automatic remeshing An effective method is presented to automatically remesh the model in numerical simulation of the fatigue crack growth in cracked panels repaired with a FRP composite patch. In the numerical simulation, both the boundary element mesh of the cracked panel and the finite element mesh of the composite patch are remeshed in each increment step of the number of cyclic loading. A typical boundary element mesh of the cross section of the cracked panel, on which the crack surface is on the left of the crack front, is shown in Fig. 3. On this cross section, only the rectangular zone a 0 bcd 0 is remeshed during the numerical simulation of the crack growth. In each increment step of the number of cyclic loading, the increment of crack length at each node on the crack front is determined based on Paris law, which is discussed in detail in Section 3.1. Using the increment of crack length and the normal direction of the current crack front at each node, the new position of each node on the crack front can be determined. The new crack front profile obtained by smoothing the curve determined with the new coordinates of the nodes on the crack front is then discretized by the nodes uniformly located along the thickness direction of the cracked panel. To ensure the accuracy of the numerical simulation, the sizes of the two column elements
Traction Singular Quarter-Point Element b
a
b’
a’
c
d
c’
d’
Crack Front
Fig. 3. Boundary element mesh of cross section with crack surface.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2555
on each side of the crack front, i.e., the elements in the zone adc 0 b 0 , are smaller than those of the others and keep the same during the simulation of the fatigue crack growth. It is noted that the coordinates of the nodes along the lines bc and a 0 d 0 , as shown in Fig. 3, keep the same during the numerical simulation of the crack growth, and the positions of the nodes along edges ad and b 0 c 0 change with the growth of the fatigue crack. In addition, the new positions of the other nodes in zones abcd and a 0 b 0 c 0 d 0 are automatically generated at each increment step of the number of cyclic loading by means of isoparametric transformation, and a bias factor is used to ensure the elements close to the crack front being smaller than those far away from the crack front. On the top and bottom surfaces of the cracked panel, the boundary elements which include the nodes along the edges ba 0 and cd 0 are modified in each step. In those elements, the new positions of the midpoint nodes of each element are located at the midpoints of the four edges connecting the four new corner nodes of the element. It is noted that the small elements around the crack tips on the top and bottom surfaces of the cracked panel keep their shapes and move with the crack tips in the numerical simulation. Meanwhile, the boundary element mesh of the other part of the cracked panel keeps the same in each step. After remeshing the boundary element portion of the model in an increment step of the number of cyclic loading, the finite element mesh of the composite patch is updated simply by replacing the x and y coordinates of the finite element nodes with the new x and y coordinates of the corresponding boundary element nodes on the attachment surface of the cracked panel. Some numerical examples have been used to demonstrate that this remeshing method is effective for this typical problem. It is noted that if the shapes of the boundary elements on the top and bottom surfaces of the cracked panel, which include the nodes along the edges ba 0 and cd 0 as shown in Fig. 3, and the corresponding finite elements of the composite patch become irregular due to a long distant propagation of the crack, a manual remeshing is needed at this step.
3. Numerical models and experiments 3.1. Numerical models The fatigue crack growth behavior of a cracked aluminum panel repaired with a circular FRP composite patch, as shown in Fig. 1, is numerically simulated. The material of the cracked panel is 2024-T3 aluminum and the adhesive is AF-163-2K. The composite patch is a ½0 3 = 45 =0 2 s glass/epoxy laminate, where the fiber direction of 0 is parallel to the longitudinal direction of the panel. The dimensions of the cracked aluminum panel, FRP composite patch and adhesive layer are listed in Table 1, and the corresponding mechanical properties are shown in Table 2. The cracked aluminum panel is subject to a cyclic stress with amplitude Dr0 = 59.6 MPa, and the applied maximum stress is rmax = 79.4 MPa and the stress ratio R = 0.25. The initial crack length on the surface of this cracked panel repaired with composite patch is 2a = 21.30 mm. In the numerical simulation, the taper of the patch, as shown in Fig. 1, is ignored and its diameter is assumed to be D. The adhesive layer is simplified as linear springs as discussed in Section 2.3. To reduce
Table 1 Dimensions of cracked panel with a FRP composite patch Components
Dimensions
Cracked aluminum panel FRP composite patch Adhesive layer
L = 248 mm, W = 118 mm, tp = 2 mm D = 50 mm, tl = 3.2 mm ta = 0.1 mm
2556
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
Table 2 Mechanical properties of cracked panel with a FRP composite patch Materials
Mechanical properties
2024-T3 aluminum Lamina of composite AF-163-2K adhesive
E = 72.39 GPa, m = 0.33 E11 = 44.12 GPa, E22 = 9.65 GPa, G12 = 4.13 GPa, m12 = 0.28 Ea = 1.18 GPa, ma = 0.34
the computational cost, an equivalent transversely isotropic single layer is employed to represent the original FRP composite patch. In order to obtain the effective material properties of this equivalent single layer, an optimization problem has been set up, in which the objective function is chosen as the sum of the leastsquares of differences between the in-plane stiffness and bending stiffness of this single layer and those of the original composite laminate, and the design variables are the material constants. Through the minimization of this objective function, the material constants of the FRP composite patch represented by the equivalent single layer are obtained as E11 ¼ 37.35 GPa;
E22 ¼ 11.38 GPa;
G12 ¼ 5.97 GPa;
m12 ¼ 0.38
where subscripts 1 and 2 are respectively parallel to the longitudinal and transverse directions of the panel. Owing to the symmetry of the problem, only one-fourth of the model is discretized. The boundary element mesh and finite element mesh are shown in Fig. 4. For the boundary element portion of the model, at the symmetric surface x = 0, the nodal displacements ux = 0, and at the symmetric surface y = 0, except for the crack surface, uy = 0. In the finite element portion, the symmetric boundary conditions are set at the two symmetric edges of the composite patch and the displacement in z direction at the top of the node at the center of the circular patch is set to zero. The average stress is applied on the loading end of the cracked panel. It is noted that the dimension in thickness direction of the schematic model in Fig. 4 is amplified by three times for clarity of view. It is important to investigate how many boundary element layers in the thickness direction of the cracked panel shown in Fig. 4, around the crack front area, are adequate to get accurate stress intensity factors along the crack front. Several mesh models of the unpached cracked panel with different numbers of boundy
FRP Composite Patch
z
Aluminum Panel
Crack Surface Crack Front
x
Fig. 4. BEM/FEM mesh of cracked aluminum panel repaired with a FRP composite patch.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2557
Table 3 Stress intensity factors of cracked aluminum panel with a FRP composite patch z (mm)
KI (MPa m1/2) with BEM/FEM
KI (MPa m1/2) with FEM
Error (%)
0.40 0.80 1.20 1.60
12.26 10.98 9.46 7.72
12.60 11.09 9.52 7.84
1.72 1.00 0.63 1.59
ary element layers in the thickness direction of the cracked panel are analyzed by means of the BEM. After comparison of these results, it is concluded that five element layers in the thickness direction of the cracked panel are sufficient to ensure the accuracy of the numerical simulation in this case. To demonstrate the validity of the combined BEM/FEM and the mesh of the model, as shown in Fig. 4, a model with a straight crack of the length a = 15 mm is analyzed by means of the combined BEM/FEM and the FEM. The mesh of the combined BEM/FEM model is similar to that as in Fig. 4. The mesh of the FEM model is much fine, which is composed of 61,046 nodes and 13,292 three-dimensional 20-node isoparametric elements, and the singular quarter-point elements are used to model the crack front area. The stress intensity factors at the points corresponding to the boundary element corner nodes along the crack front, except for the cross points of the crack front with the top and bottom surfaces of the cracked panel, obtained with these two methods are compared in Table 3. It is shown that the difference between the stress intensity factors obtained with the combined BEM/FEM and the FEM are small. The agreement of the stress intensity factors at the points on the crack front demonstrates the validity of the combined BEM/ FEM and the mesh model used in the numerical simulation of fatigue crack growth in the cracked aluminum panel repaired with a FRP composite patch. In the numerical simulation of the fatigue crack growth, the stress intensity factor KI max of the cracked panel with a certain crack front profile under the action of maximum stress rmax is firstly determined by means of the combined BEM/FEM. Because the stress intensity factors at the cross points of the crack front with the top and bottom surfaces of the cracked panel cannot be obtained accurately, the stress intensity factors at these points are modified by smoothing with the values at other nodes along the crack front. Taking account of Eqs. (3) and (4), the stress intensity factor range DKI corresponding to the stress range Dr0 can be determined by DK I ¼ ð1 RÞK I max
ð16Þ
The increment of crack length at each node along crack front, corresponding to an increment of the number of cyclic loading, can then be determined based on Paris law da ¼ CðDK I Þm dN
ð17Þ
where da/dN is the crack growth rate, and C and m are the material constants of the aluminum panels. In addition, because the influence of residual stress generated in the process of patching FRP composite laminates on stress intensity factor range is not great [1], the residual stress is neglected in the numerical simulation. The disbond effect of patch during fatigue crack growth process is also ignored in the numerical simulation for simplification. 3.2. Experiments To determine the material constants of fatigue crack growth of 2024-T3 aluminum, fatigue test of an unpatched center-cracked aluminum panel was performed. In this test, the cracked aluminum panel was subject to a sinusoidal stress with amplitude Dr0, a stress ratio of 0.25 and a cyclic frequency of 0.5 Hz
2558
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
Table 4 Dimensions and loading of unpatched cracked aluminum panel specimen L (mm)
W (mm)
tp (mm)
2a (mm)
Dr0 (MPa)
248
118
5.0
60.0
23.8
40
a (mm)
35
30
Numerical result Experimental data
25
20 0
1
2
3
4
[×104]
N (cycles) Fig. 5. Crack length of unpatched cracked aluminum panel versus number of cyclic loading.
in the environmental temperature of 27 C. The dimensions of the specimen, the initial crack length 2a and the applied cyclic stress Dr0 are listed in Table 4. The experimental data of the crack length versus the number of cyclic loading are shown in Fig. 5. Fatigue test of the cracked aluminum panel repaired with a circular FRP composite patch, as shown in Fig. 1, was then performed. The materials and dimensions of the cracked panel, adhesive and composite patch, as well as the applied stress are the same as those of the model discussed in Section 3.1. The other test conditions are the same as those in the test of the unpatched cracked aluminum panel mentioned above. Before the cracked aluminum panel is repaired, it was applied some cycles of loading until the crack length on the surface arrives at 2a = 21.30 mm. The cracked aluminum panel was then repaired with a circular FRP composite patch and was continuously applied the same cyclic loading until fatigue failure. The experimental data of the crack length on the unpatched side of the cracked panel repaired with a FRP composite patch versus the number of cyclic loading are shown in Fig. 6. For the purpose to investigate the fatigue crack growth behavior, the specimens of the unpatched and repaired aluminum cracked panels were fractured after some cycles of loading, and the crack front profiles of the specimens were then observed. The micrographs of typical examples of the crack front profiles corresponding to these two cases are respectively shown in Figs. 7 and 8.
4. Numerical results and discussion Fitting the results of the numerical simulations of fatigue crack growth of the unpatched cracked panel with the BEM to the experiment data gives the material constants of 2024-T3 aluminum C = 7.04 · 107
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2559
25
∆σ0=80.0 MPa 20
∆σ0=59.6 MPa
a (mm)
15
∆σ0=40.0 MPa
10
Present numerical model Simplified numerical model Experimental data
5
0 0
1
2
3
4
5
6
7
8
9
10 [×104]
N (cycles) Fig. 6. Crack length on unpatched side of cracked aluminum panel versus number of cyclic loading.
Fig. 7. Crack front of unpatched specimen (material: 2024-T3 aluminum; thickness: 5 mm; number of cyclic loading: 46448).
(MPa m1/2)2.48 mm/cycle and m = 2.48. In the numerical simulations, the mesh of the boundary element model is similar with those of the boundary element portion of the combined BEM/FEM model as shown in Fig. 4. Besides the boundary conditions as set in the boundary element portion of the model in Fig. 4, the displacement in z direction of the node located at the origin of the coordinate system is set to zero in this BEM model. The variation of crack length a versus the number of cyclic loading N determined by the numerical simulation is illustrated in Fig. 5, and it agrees well with the experimental data. In addition, the crack front profiles and the distributions of stress intensity factors of the unpatched specimen under the maximum stress rmax = 31.76 MPa at some cycles of loading during the fatigue crack growth simulated with the BEM are shown in Fig. 9. It is obvious that with the change of the distribution of stress intensity factors along crack front, the crack front is gradually curved from the initial straight line and propagates nearly parallel after some cycles of loading. The shapes of the crack fronts in the process of stable propagation determined by numerical simulation as in Fig. 9(a) is similar to that observed in experiment as in Fig. 7.
2560
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
z (mm)
Fig. 8. Crack front of specimen repaired with a FRP composite patch (material: 2024-T3 aluminum; thickness: 2 mm; number of cyclic loading: 61,600).
5 4 3 2 1 0 22
A
24
26
E
D
C
B
28
30
32
F
34
36
G
38
40
42
x (mm)
(a) 20.0
G F
15.0
KI (MPa m1/2)
E D 10.0
A
5.0
B
C
0.0 0.0
(b)
1.0
2.0
3.0
4.0
5.0
z (mm)
Fig. 9. Crack fronts and stress intensity factors in the case of unpatched cracked aluminum panel (number of cyclic loading: A: 0; B: 1000; C: 8000; D: 25,000; E: 49,000; F: 66,000; G: 78,000). (a) Crack fronts and (b) stress intensity factors.
Using the material constants, the combined BEM/FEM is employed to numerically simulate the fatigue crack growth behavior of the cracked aluminum panel repaired with a FRP composite patch. It is noted that the fatigue crack growth process of the cracked aluminum panel before being repaired is analyzed with
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2561
the BEM, and the final crack front profile in this BEM simulation is adopted as the initial profile of the repaired cracked panel. The variation of crack length on the unpatched side of the cracked aluminum panel repaired with a FRP composite patch, in the case of Dr0 = 59.6 MPa, versus the number of cyclic loading, which is determined by the numerical simulation, is demonstrated in Fig. 6. It is shown that the curve determined with the numerical method agrees well with the experiment data. The numerically simulated crack front profiles and distributions of stress intensity factors of the cracked aluminum panel under the peak stress rmax = 79.4 MPa are illustrated in Fig. 10. After the cracked aluminum panel is repaired, the deformation of the cracked panel is constrained by the FRP composite patch, so that the stress intensity factors along crack front decrease. The constraint of the FRP composite patch is implemented through the shear force transferred by adhesive layer. At the very beginning of fatigue crack growth after its reinitiating, the cross point of crack front with the unpatched surface of the cracked panel grows quickly due to the larger stress intensity factor, comparing with other points along the crack front. Meanwhile, the stress intensity factor at the cross point decreases and those at the other points increase as the number of cyclic loading increases. Therefore, the growth rate of crack on the unpatched side decreases slightly during the initial cyclic loading. By comparing Figs. 8 and 10(a), it is found that the crack front profiles determined with the combined BEM/FEM numerical simulation are similar to that observed in experiment except the portion close to the unpatched side of the cracked panel. As aforementioned, most of the previous works such as Ref. [12] did not consider the variation of stress intensity factors along curved crack front. However, as illustrated in Figs. 8 and 10, the crack front profiles of
z (mm)
2.0 1.5
B
1.0
A
0.5 0.0 10
D
C 11
12
13
14
F
E 15
16
17
G
18
19
20
x (mm)
(a) 15.0
F
10.0
E
D
1/2
KI (MPa m )
G
C B
5.0
A
0.0 0.0
(b)
0.5
1.0
1.5
2.0
z (mm)
Fig. 10. Crack fronts and stress intensity factors of cracked aluminum panel repaired with a FRP composite patch (number of cyclic loading: A: 0; B: 4000; C: 8000; D: 30,000; E: 45,000; F: 66,000; G: 82,000). (a) Crack fronts and (b) stress intensity factors.
2562
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
cracked panel repaired with a FRP composite patch are not straight lines and the distributions of the stress intensity factors along the crack front are nonlinear during the fatigue crack growth. Therefore, it is important to consider these two factors in the numerical simulation of the fatigue crack growth behavior of the cracked panel repaired with a FRP composite patch. Using the average of stress intensity factors along the crack fronts, except the values at cross points of the crack front with the top and bottom surfaces of the cracked panel, and assuming that the crack front keep straight and perpendicular to the top and bottom surfaces of the cracked panel and propagate parallel during the crack growth, the fatigue crack growth of the same model as shown in Fig. 4, under the cyclic stress Dr0 = 59.6 MPa, is numerically simulated by means of the combined BEM/FEM. For the convenience of description, this model is referred as the simplified numerical model. The curve of the crack length versus the number of cyclic loading of the simplified numerical model is shown in Fig. 6 with a dash line. It is obvious that the curve of the simplified numerical model is close to the experimental data at the initial stage, but the error becomes larger and larger with the increase of the number of cyclic loading. Although the results of the simplified numerical model are closer to the experimental data, in the initial stage, than those of the present numerical model in which the curved crack front and the nonlinear distribution of the stress intensity factors along the crack front are taken into account, the results of the later are closer to the experimental data in the whole stage of the fatigue crack growth. To investigate the fatigue crack growth behavior of the cracked panel repaired with a FRP composite patch under two different cyclic stresses with amplitudes Dr0 = 40.0 MPa and 80.0 MPa, the numerical simulations are carried out with the combined BEM/FEM. The variations of the crack length on the unpatched side of the cracked aluminum panel versus the number of cyclic loading in these two cases are shown in Fig. 6. It is illustrated that as the increase of the amplitude of the cyclic loading, the crack front propagates much faster.
5. Conclusions Fatigue crack growth behavior of a cracked aluminum panel repaired with a FRP composite patch under uniaxial cyclic loading is numerically investigated by the combined BEM/FEM. The curve of fatigue crack length on the unpatched side of the cracked panel versus the number of cyclic loading, which is determined by the numerical simulation, agrees well with experimental data. The crack front profiles and the distributions of stress intensity factors along the crack fronts during fatigue crack growth are numerically simulated. All these results demonstrate that the numerical method used in this paper is efficient in the evaluation of fatigue crack growth behavior of cracked aluminum panels repaired with an adhesively bounded FRP composite patch.
Acknowledgements This work was partly supported by the Grand-in-Aid for Scientific Research (A) (2) No. 14205138 and (B) (2) No. 13555264 from the Ministry of Education, Culture, Sports, Science and Technology of Japan to HS. The authors would like to thank Mr. Takamitsu Sato for his assistance to conducting the experiments.
References [1] Baker AA. Repair efficiency in fatigue-cracked aluminum components. Fatigue Fract Engng Mater Struct 1993;16(7):753–65. [2] Baker AA, Rose LRF, Jones R. Advances in the bonded composite repairs of metallic aircraft structures, vols. 1&2. London: Elsevier Science Publications; 2002.
H. Sekine et al. / Engineering Fracture Mechanics 72 (2005) 2549–2563
2563
[3] Denney JJ, Mall S. Characteristics of disbanded effects on fatigue crack growth behavior in aluminum plate with bonded composite patch. Engng Fract Mech 1997;57(5):507–25. [4] Schubbe JJ, Mall S. Investigation of a cracked thick aluminum panel repaired with a bonded composite patch. Engng Fract Mech 1999;63(3):305–23. [5] Seo DC, Lee JJ. Fatigue crack growth behavior of cracked aluminum plate repaired with composite patch. Compos Struct 2002;57(1–4):323–30. [6] Ratwani MM. Analysis of cracked adhesively bonded laminated structures. AIAA J 1979;17(9):988–94. [7] Sun CT, Klug J, Arendt C. Analysis of cracked aluminum plates repaired with bonded composite panels. AIAA J 1996;34(2):369–74. [8] Naboulsi S, Mall S. Modeling of a cracked metallic structure with bonded composite patch using the three-layer technique. Compos Struct 1996;35(3):295–308. [9] Naboulis S, Mall S. Nonlinear analysis of bonded composite patch repair of cracked aluminum panels. Compos Struct 1998; 41(3–4):303–13. [10] Schubbe JJ, Mall S. Modeling of cracked thick metallic structure with bonded composite patch repair using three-layer technique. Compos Struct 1999;45(3):185–93. [11] Naboulsi S, Mall S. Fatigue crack growth analysis of adhesively repaired panel using perfectly and imperfectly composite patches. Theor Appl Fract Mech 1997;28(1):13–28. [12] Shibuya Y, Fujimoto S, Aoki D, Sato M, Shirahata H, Fukunaga H, et al. Evaluation of crack growth in cracked aluminum panels repaired with a bonded composite patch under cyclic loading. Adv Compos Mater 2001;10(4):287–97. [13] Kam TY, Chu KH, Tsai YC. Fatigue of cracked plates repaired with single-sided composite patch. AIAA J 1998;36(4):645–50. [14] Chue CH, Chang LC, Tsai JS. Bonded repair of a plate with inclined central crack under biaxial loading. Compos Struct 1994;28(1):39–45. [15] Young A. Three-dimensional analysis of a patched cracked solid using a combined boundary element and finite element model. DRP 26C-R02-Milestone 2 (Item 410), Working Paper MS4-93-WP-24, Materials and Structures Department, Defense Research Agency, Hampshire, 1993. [16] Brebbia CA. The boundary element method for engineering. London: Pentech Press; 1978. [17] Pagano NJ. Stress field in composite laminates. Int J Solids Struct 1978;14(5):385–400.