Crack Propagation in Fiber-Reinforced Composites

(1)

TVE 13 032 juni

Examensarbete 15 hp

Juni 2013

Crack Propagation in Fiber-Reinforced

Composites

A semi-analytical method for determining

cohesive laws

Linda Storvall

David Bergstrand

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Crack Propagation in Fiber-Reinforced Composites - A

semi-analytical method for determining cohesive laws

Linda Storvall, David Bergstrand

When describing material properties in fiber-reinforced composites, cohesive laws are needed. Cohesive laws can be thought of as a nonlinear spring trying to

counteract the crack propagation. In the laboratory two different methods were used to gather data to determine these laws. One with a complex rig, applying pure moments to the test specimen, but with a simple mathematical analysis and one with a simple rig, applying pure tension to the specimen, but with a complex mathematical analysis. To approach the cohesive laws in the pure moment tests only experimental data was needed. In the tension test a semi-analytical approach, using both

experimental data and purely theoretical formulas, was required. The data and formulas was implemented in MATLAB. The computations for the pure tension tests became so time consuming that no results were calculated within the given time frame, even when a supercomputer was used.

Handledare: Kristofer Gamstedt

(3)

Populärvetenskaplig sammanfattning

En komposit är ett material där era olika typer av material blandas för att nå önskade egenskaper.

Det kan vara allt från att förstärka betong med gjutjärn för att få en starkare konstruktion till att blanda olika typer av brer med plast, så som kolber, för att nå önskade egenskaper. Kompositer där brer och plast blandas kallas ofta för berkompositer, vilka ses lite som ett framtidens material och används i t.ex ygplan och rymdfarkoster. För att ta reda på ett materials egenskaper måste det testas. I det här projektet testades två olika metoder för att mäta hur starkt en typ av berkomposit är. Det ena testet är bra på så sätt att matematiken för att komma fram till materialets egenskaper är lätt, medan testet innebar en svår konstruktion. Det andra testet är helt tvärtom, lätt att genomföra, men kräver mycket beräkningar för att komma fram till egenskaperna. Det visade sig vara så tunga beräkningar att inga resultat uppnåddes, inte ens när beräkningarna utfördes på en superdator.

(4)

1 Introduction

The study is based on an unpublished article by E.K. Gamstedt and S.P. Fernberg [1]. All the assumptions and statements referenced to this article are directly taken from the article by consultation with the authors.

1.1 Background

A composite material is made out of two or more materials with dierent physical and/or chemical properties. This results in a material in which the properties can dier a lot due to the individual components in it. Dierent kinds of composites are for example concrete, paper and ber-reinforced composites such as berglass, metal composites and ceramics composites [5]. In this project only

ber-reinforced composites are concerned. The specimens in the tests are made of berglass [1].

Composites can be composed in a vast amount of ways to receive wanted properties of the material.

In cases where the material has to carry a weight or in other ways deal with external forces, a high fracture toughness is desired. A ber-reinforced composite has a very high fracture toughness due to pull-out of bers just behind the crack tip. Fiber pull-out is when the bers of the material, instead of failing at the crack tip, are being pulled out as bridges over the crack, as shown in g. 1.1, and eventually fail because of elongation of the bers. The problem with ber-reinforced composites is that the inelastic damage zone at the crack tip is large and that makes classical elastic fracture mechanics a poor model for describing crack propagation. Due to these special characteristics of ber- reinforced composites the crack propagation is more accurately described with cohesive laws. When pull-out occurs, the bers being pulled out create friction against the composite material, and therefore functions as a resistance to the crack propagation. This resisting force, the closing traction, can be modelled as nonlinear springs holding the crack together just behind the crack tip, as shown in g.

1.1.

Figure 1.1: A comparison between the physical cracking and the cohesive model where the bridging

bers work as springs.

The cohesive law is expressed as the closing traction between opposite crack surfaces, σ, and is a function of the crack opening displacement, δ. The crack opening displacement can be explained as the opening at the apex of the crack, or as in this case, the opening just behind the pulled out

bres. In the model using nonlinear springs, the crack opening displacement is the opening at the spring positioned furthest from the crack tip. The crack opening displacement usually depends on the position x along the crack plane, i.e. δ (x). In this study the cohesive law is assumed to be

(6)

σ (δ) =

(σ0− kδ, 0 ≤ δ ≤ δ0

0, δ ≥ δ0

(1.1) where σ0is the maximum closing traction andδ0is the critical crack opening displacement, i.e. the maximum crack opening just before the bers break. The constant k is expressed as

k = 2 (σ₀δ₀− Γ)

δ²₀ (1.2)

[1], where Γ is the fracture energy.

This study centers around determining the parameters of the cohesive law, using a semi-analytical method, and making a comparison between dierent methods in determining these parameters. The semi-analytical method is a combination of using experimental data and theoretical formulas. In the following chapters, the experiments and analyses are presented.

1.2 Laboratory experiments

Prior to this study, the cohesive laws were experimentally measured with two dierent methods. Both methods used a double-cantilever beam (DCB) for testing. The advantage with a DCB-test is that it has a simple geometry, as presented in g. 1.2.1, which makes the testing easy and straightforward.

One method for measuring the cohesive laws in the laboratory was by loading the DCB-test with pure moments and the other method was by loading the test with wedge forces, creating pure tension.

Although a DCB-test has a simple geometry, the crack can propagate chaotically. To get a good measurement the crack has to be controlled. This was done by making a small groove through the middle of the specimen where the crack propagation is wanted. The specimen was made out of a sheet moulding compound (SMC) based on an unsaturated polyester resin with composition according to Table 1 (Bopreg 104, Reichhold). The material is reinforced with glass bres [1].

Figure 1.2.1: Geometry och the double-cantilever beam.

The procedure of the two dierent kinds of experiments dier somewhat in diculty. A pure tension DCB test consists of applied wedge forces, pulling the edges of the specimen straight upwards respectively straight downwards, creating a crack propagation. The rig is presented in g. 1.2.2. This experiment can easily be done in any solid mechanics laboratory with access to a tensile test machine.

Pure moment tests, however, require a more complex structure. The apparatus has to be modied with a rig transforming pure tension from the machine to pure moments on the specimen. This rig has to be carefully designed not to accidently apply tension to the specimen. The structure of the moment loaded experiment can be seen in g. 1.2.3. Due to this dierence in diculty, pure tension tests are

(7)

preferable to execute. However, there is also a considerable dierence in describing the two dierent methods analytically.

Figure 1.2.2: Photograph of the pure tension test rig.

Figure 1.2.3: Photograph of the pure moment test rig.

(8)

2 Theory

2.1 The J-integral approach

Using the J-integral is one approach to describe crack propagation in materials. This integral denes the energy release rate with respect to the surface area [6]. By plotting the J-integral against the crack root opening displacement, δ^∗, a description of the crack propagation is given. This method accounts for the cohesive laws directly, since the cohesive law can be expressed as the derivative of the J-integral with respect to the crack root opening displacement [7, 8]. In g. 2.1 the formulas for calculating the J-integral for pure moments and pure tension is presented.

Figure 2.1: Picture of the dierent DCB testing methods and their solution to the J-integral.

On rst sight this method might seem easy. However, the angle θ is required for calculating the J-integral for wedge loaded test. θ is very dicult to measure in a laboratory experiment, and is therefore seldom gathered as data. For testing with pure moments the solution of the J-integral is a function depending only on the geometry of the specimen and the applied moment, M, as

J = 12(1 − ν²)M²

BbH³E , (2.1)

where B is the beam width, b is the notch width and H is the half height of the beam. The J-integral is therefore a very straightforward method of calculating the cohesive law for pure moment

(9)

tests. For pure tension tests an alternative method is required. In this case a semi-analytical method was tested.

2.2 The R-curve approach

In materials where the damage zone is small compared to the dimensions or the structure of the crack, the crack propagation is easily described by the relation between the energy release rate and the crack length [9]. This is called an R-curve, as seen in g. 2.2. However in materials, such as ber-reinforced composites, the damage zone is not small compared to the dimensions or the structure of the crack.

This means, in addition to the R-curve, a cohesive law has to be determined [1]. Since the R-curves depend on the geometry of the material it is not a characteristic of the material, on the other hand the cohesive laws depend only on the material characteristics and is therefore a way to describe material properties. To nd the material properties the cohesive law is usually described by the cohesive crack surface traction with respect to the crack opening displacement, σv(δ) [1].

Figure 2.2: Example of a tted simulated R-curve from experimental data.

For tests concerning wedge forces the R-curve approach is a better method to attack the problem [1]. The R-curve from a DCB test depends on the measured crack length, ∆a, and the measured wedge force, P . This means that the R-curve needs to be discretized for every measured crack length and wedge force. An initial crack is made to the specimen and is denoted a0.

GR(∆ai, Pi) =2Pi(a0+ ∆ai)²

bE⁰I . (2.2)

From eq. 2.2 the R-curve is easily shown by plotting GR with respect to the crack length ∆a.

When calculating a theoretical R-curve the cohesive laws aren't considered and a dierent approach is therefore needed. The stress intensity approach is used to give an accurate R-curve that takes the cohesive laws in consideration [10]. The stress intensity factor can be written as

K₀= K_a(∆a) − K_b(∆a), (2.3)

where Ka is the stress intensity from the wedge forces, Kb is from the cohesive action and K0 is the stress intensity factor. Since the cohesive traction only works just behind the crack tip a new

(10)

variable, x, is dened to denote the cohesive action. This variable goes along the crack length and is therefore dened by 0 < x < ∆a. This means the cohesive contribution is dened as σ(x). However, to let the cohesive traction account for R-curve behaviour and let the cohesive laws work for more than one geometry, the cohesive laws are dened by the local crack opening σ(δ(x, ∆a)). To denote the cohesive contribution in the stress intensity approach a weight function, m(δ(x, ∆a)), is needed.

Since the cohesive traction varies along the crack length an integral is needed to sum up all the small cohesive actions [1]. From these concerns a new relation can be obtained and is written as

δ(x, ∆a) = δ_a(x, ∆a) − 1 E⁰

ˆ ∆a 0

σ(δ(ξ, ∆a))m^∗(x, ξ, ∆a)dξ, (2.4) where E' is dened by

E⁰= E

1 − ν², (2.5)

for plane strain, where E and ν are the isotropic Young's modulus and Poisson ratio, respectively.

The weight function for δ(x, ∆a) can be determined to be

m^∗(x, ξ, ∆a) = ˆ ∆a

ζ

m(ξ, α)dα, (2.6)

where

ζ = max(x, ξ) (2.7)

[14].

Eq. 2.4 is very similar to a Fredholm integral equation of the second kind and has to be solved numerically. The only dierence to a normal Fredholm equation is that the crack opening displacement is sought and is also an argument in the cohesive law operator in the integral. This nonlinear equation can be solved by discretizing the crack opening displacement and then using an iterative method to

nd a function for δ(x, ∆a). Since it's so similar to a Fredholm equation of second kind it is guaranteed to converge as long as the cohesive zone is much smaller than the crack length [13].

From the weight function, the crack opening displacement and a reference stress intensity factor the nominal crack opening displacement can be expressed by

δa(x, ∆a) = P E⁰

ˆ ∆a x

m(x, α)m(−a0, α)dα. (2.8)

The weight function can be expressed by a number of ways. A fast and straightforward weight function can be determined from Bernoulli-Euler beam theory. The weight function then becomes

m₁(x, ∆a) =

√12(∆a − x)P

H^3/2 (2.9)

[2] However, the accuracy of eq. 2.9 might be questionable and a more precise weight function might be to prefer. Fett and Munz determined a more precise weight function by a boundary collocation method and is written as

m₂(x, ∆a) = r12

H

∆a − x H + 0.68

+

s 2

π(∆a − x)exp (

−

r12(∆a − x) H

)

(2.10)

(11)

[3].From a two-dimensional elasticity and a Wiener-Hopf technique a third weight function [4] with even higher accuracy can be achieved

m3(x, ∆a) = r12

H

∆a − x

H + 0.673

+

s 2

π(∆a − x)− 1 H

(

0.815 (∆a − x) H

⁰

.619 + 0.429 )−1

. (2.11) From (Ksumma) the stress intensity factor can be rewritten as a function of only the intensity of the wedge forces

Ka(∆a) = K0+ ˆ ∆a

0

m(x, ∆a)σ(δ(x, ∆a))dx (2.12)

which is related to the applied energy release rate as [15]

Ga(∆a) = K_a(∆a)²

E⁰ . (2.13)

The relation between the measured stress energy release rate, GR, and the applied stress energy release rate, Ga, makes it possible to t a regularized cohesive law [1]. The tting is done by

χ²(c, J0) = min

c,J₀ n

X

i=1

{Ga(∆a_i, c, J₀) − G_R(∆a_i, P_i)}²

n , (2.14)

where c is a vector containing σ0, δ0and Γ from eq. 1.1. The minimization of eq. 2.14 with respect to σ0, δ0, Γ and J0 gives the parameters that are sought in the cohesive law. However, J0 can be estimated as J0= G_R(∆a = 0)and therefore only σ0, δ0and Γ needs to be found in the minimization.

The stress intensity factor K0is related to J0 by

K₀²= J₀E⁰ (2.15)

By using this method to nd the cohesive law, it's easy to change material, geometry or forces in testing and still nd the parameters of the cohesive law. As long as an R-curve can be measured, the cohesive laws can be found [1].

3 Methods

To nd the parameters of the cohesive laws, the minimization of eq. 2.14 has to be found by comparing the theoretical energy release rate with the experimental energy release rate from the wedge loaded tests. This is done by separating the dierent formulas to a number of functions in MATLAB. When the sought after parameters are found, the resulting cohesive law can be compared to the experimental results from the moment loaded tests.

The data given from the pure moment experiments is a set of crack root opening displacements, δ^∗, with its respective moments, M. For the pure tension tests the wedge forces, P , and the crack length,

∆a, was given. Seven dierent moment loaded DCB tests and two dierent wedge loaded DCB tests has been performed on the same type of material.

(12)

3.1 Implementing the J-integral approach

Using eq. 2.1, a graph of the J-integral's variety over the crack root opening displacement is made from the pure moment experimental data. Approximating this experimental J-integral curve to a polynomial and dierentiating it with respect to the crack root opening displacement will result in the experimentally given cohesive law for the pure moment tests.

3.2 Implementing the R-curve approach

By implementing the R-curve approach for the wedge loaded tests, the cohesive law is eventually found.

To begin with, the given data for the wedge loaded tests, ∆a and P , are loaded in as input parameters.

From this data the experimental energy release rate, GR, can be found using eq. 2.2. An R-curve is given by plotting a graph of the experimental energy release rate with respect to the crack length.

This R-curve is tted to a low-order polynomial to nd a good approximation for J0, at GR(∆a = 0). By doing this, J0is no longer an unknown parameter in eq. 2.14. The remaining unknown parameters are nowσ0, δ0 and Γ, all from eq. 1.1 and eq. 1.2. The minimization in eq. 2.14 is done by letting the computational program in MATLAB try out values of the parameters, near and around a well-chosen guess, and nding the parameter values that give the minimum value of χ². Guessing many dierent parameters might get unstable, therefore the assumption

k = σ0

2δ0 (3.1)

is made, to replace eq. 1.2. Now the unknown parameter Γ is eliminated from eq. 1.2 and only two unknown parameters remain, σ0 and δ0. To get a reasonable guess, σ0 is set to the maximum shear stress of the material, σmax, and δ0is determined by a visual estimate from the experiment. The method was implemented for all three dierent weight functions, eq.2.9, eq. 2.10 and eq. 2.11.

3.3 Comparing the methods

When the estimated parameters are found, the cohesive law is determined for the pure tension tests and is set in comparison to the cohesive law from the pure moment tests. This is either done by dierentiating the J-integral from the pure moment tests, as earlier mentioned, or by integrating the cohesive law from the pure tension tests and comparing the result to the J-integral. The latter is a more stable method, since dierentiating might enhance noise from the data.

4 Results

4.1 J-integrals

The J-integrals with respect to the crack length from the seven dierent pure moment DCB tests were plotted in the same graph, as in g. 4.1. They purely depend on the measured data from the experiments. The J-integral's dependence on the data is presented in eq. 2.1.

(13)

−1 0 1 2 3 4 5 0

0.5 1 1.5 2 2.5 3 3.5x 10⁶

Crack root opening displacement [m]

J−integral [J/m2]

Figure 4.1: Graph showing the J-integral with respect to the crack root opening displacement from the seven pure moment tests.

4.2 R-curves

From the two pure tension tests one R-curve each, depending solely on the data from the experiments, was plotted, as presented in g. 4.2. These are modelling the energy release rate in eq. 2.2 with respect to the crack length. Since there was an initial crack at a0 = 0.1mm from the edge made in the untouched specimen, before the tests were performed, the crack length does not begin at zero.

A polynomial curve, approximated to the experimental R-curves, was plotted alongside the R-curves.

Where these polynomial curves cross the initial crack length, the estimated value for J0can be found.

(14)

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.8

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8x 10⁴

Crack length [m]

Stress energy release rate [J/m2]

0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.8

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8x 10⁴

Crack length [m]

Stress energy release rate [J/m2]

Figure 4.2: Graphs showing the R-curves gathered from data with their tted curves. The tted curves are 4th grade polynomials.

4.3 Cohesive law for the pure tension tests

The computations were very time consuming and therefore no results were given for the parameters σ₀ and δ0within the given time frame. Consequently, no cohesive law was found for the pure tension tests.

4.4 Comparison

The J-integral from the pure moment tests could not be compared to the resulting J-integral from the pure tension tests, since the cohesive law was not found for the pure tension tests, due to time consuming computations.

5 Discussion

5.1 Choosing a method for integrating in MATLAB

When we implemented the equations in MATLAB, the chosen method for calculating integrals was MATLAB's built-in method integral. This method calculates the integrals to a high order of accuracy (to the tenth decimal) and is MATLAB's latest built-in method for calculating integrals [11]. However, this method wasn't the only one tested.

To calculate an integral with a high order of accuracy without too many steps we tried using Simpson's method. Simpson's method was implemented in two ways, both with MATLAB's built-in quad and a integral solver written by us. We soon found that quad calculated much faster than our own code, and quickly abandoned the idea of writing our own integral solver. quad was also somewhat faster than integral, but we found that quad had problems handling singularities and fast convergences[12], which became a problem when we tested the program for the weight functions m2and m3.

(15)

The built-in method int was also tested. This method handles the variable of integration sym- bolically, which resulted in fast calculations. However, we found that not all of these variables could be treated as symbols, since x and ξ also functioned as boundaries to integrals, which has to take numerical values. Ultimately we chose integral for calculating the integrals numerically, as this was the method with highest accuracy and the only one able to perform our nested integral calculations correctly.

5.2 Thoughts on why the cohesive law was not found

As shown in chapter 2.2, the function χ² from eq. 2.14 ultimately consisted of both double and triple integrals, where one of the equations even was a integral equation (eq. 2.4). When we implemented this in MATLAB the result was a nest of integral functions. The Fredholm integral equation had to be calculated by iterations and some of the functions consisted of for loops to split up vectors to single digits, since the boundaries of the integrals could not be vectors. Consequently, some of the functions - constisting of integrations - were called for many thousand times during just one guess of the parameters in χ². We believe this is the main cause of the program running slowly.

Running the program on our personal laptops, which was the idea at rst, with enough calculation to be sure that a high accuracy was achieved, the program would have to run for many years before any results would be found. We limited the amount of iterations and guesses of parameters so the program would be nished running after 10,000 hours on our peronal laptops, which was the requirement for getting access to running the computations through a supercomputer. The supercomputer computed the program for ten days on 16 parallel processor cores, but that time frame turned out being too short.

No results were found. Since MATLAB itself is able to split up the work of running the program on several parallel processors, the result surprised us. We believe that in this case, MATLAB was not able to split up the work on several cores because the built-in method we used to guess values of the parameters and minimize χ², fminsearch, was not able to run its computations parallelly. We do know that the program worked, though, from running it with such small amount of iterations and guesses of parameters as possible. This did not give us a reliable result, however.

For further research on the subject, another method for minimizing would be recommended.

5.3 Suggestions for decreasing the running time

MATLAB might not have been the ultimate program for the task. MATLAB is an easy tool for doing advanced calculations, because of its built-in methods for many dierent types of mathematical calculations and because of how vectors and matrices are handled. However, MATLAB is not optimized for doing fast calculations and e.g. C++ or C would probably have given us the results much faster. On the other hand, C or C++ would have taken us a long time to learn and write the code in, since they do not come with as many built-in computational funstions as MATLAB. Since our code consists of nested integrals, both vector-valued and regular ones, and an iteration inside of the nest, it will always take a long time doing the calculations in MATLAB. The way we see it, not much more can be done to optimize the program's running time in MATLAB. The only thing we could do to make the code run faster was to change the number of iterations and the number of guesses in the minimization, but this would have a negative eect on the accuracy of the results. Antoher advantage of using MATLAB is that MATLAB itself is able to split up the work on several cores while running on a computer with multiple parallel processors, and no knowledge of parallel process programming is required. All in all, MATLAB is an easy tool to use, but for future implementation of the semi-analytical method we would recommend using a dierent programming tool for faster results.

(16)

6 Conslusions

The semi-analytic method for determining the cohesive law from a pure tension DCB test is possible to implement in MATLAB, although the computations are very time consuming. The computations were found to be so time consuming that no results were found within the given time frame. MATLAB is an easy tool for implementing advanced mathematical calculations like these, but we would recommend using a programming tool like C++ or C if time is considered a more important factor than the simplicity of writing the program.

(17)

References

[1] A Comparison of Cohesive Laws Determined from Double-Cantilever Beam Tests with Wedge Forces and Pure Moments by E.K. Gamstedt and S.P. Fernberg

[2] Broek D. Elementary engineering fracture mechanics. Dordrech: Martinus Nijho, 1986.

[3] Fett T. Weight function for double cantilever beam specimen. Int J Fracture 1996;80(1):R3-R8.

[4] Foote RML, Buchwald VT. An exact solution for the stress intensity factor for a double cantilever beam. Int J Fracture 1985;29(3):125-134.

[5] http://www.science.org.au/nova/059/059key.html, [Date of access: 2010-05-27]

[6] J. R. Rice, "A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks", Journal of Applied Mechanics, 35, 1968, pp. 379-386.

[7] Olsson P, Stigh U. On the determination of the constitutive properties of thin interphase layers an exact inverse solution. Int J Fract 1989;41(4):R71-R76.

[8] Suo Z, Bao G, Fan B. Delamination R-curve phenomena due to damage. J Mech Phys Solids 1992;40(1):1-16.

[9] http://www.sut.ac.th/engineering/metal/pdf/MechMet/11_Fracture%20mechanics.pdf, [Date of access: 2010-05-27]

[10] Fett T, Munz D, Thun G, Bahr H-A. Evaluation of bridging parameters in aluminas from R-curves by use of the fracture mechanical weight function. J Am Ceram Soc 1995;78(4):949-951.

[11] http://www.mathworks.se/help/matlab/ref/integral.html, [Date of access: 2010-05-27]

[12] http://www.mathworks.se/help/matlab/ref/quad.html, [Date of access: 2010-05-27]

[13] Olsson M, Giannakopoulos. R-curve model of certain ceramic composites. Fatigue Fract Eng Mater Struct 1993;16(5):539-554.

[14] Rice, JR. Some remarks on elastic crack-tip stress elds. Int J Solids Struct, 1972:8(6):751-758.

[15] http://enpub.fulton.asu.edu/concrete/ACBM_Faculty_1999/sld053.htm, [Date of access: 2010- 05-27]