... ... .. . . 0%. 020‘? 11 00 v.0l.. 001000- 00“. . ...... 0. . I . 4 I hé’. 0.0.“ .00 ...“ $1.”... 0“ 000' .13.“.0‘9“ 2. .rf00 0.... 00.20.00.000 V a ...}... «a. ‘5...) ..L.‘ .ZJ. ...-:03... a .0. 0130.1 €..~J..‘...I...0..D\.....:. 00 0. “$0000. ‘.0 0.0.20 .X.0.000.1. .03 0a a. .00 10.000A‘V0A 0 10 0 .0.“ .000370310‘0I . 5.0 .. 03. ”.00... “2.0.1.0.. . I...“ . ....00' . .30, 0 39"”..01‘0 lip-'0”? 30:30.30... 300.0.‘105000001' 000000910 33.0.0... .‘0. 30 0 U ’00WA0 \0 00” ‘0... 0. 000 1.000. 00.000‘1 0a.000400..000‘0..10.030000000000. V0.0..C’.( . .3000. v” .0, 0000. .0..’.\0500000 (0.020 0. ‘00. 00‘... 00... '“0.‘.0000'0O00\:. O '1! ..3000 0000030 92 .‘5b 0‘30. 08..J3..:. ... 34.00.103.000: .130! 0.00..0.0G..0000C ..lw 1‘ 3... 0.00.01.33 ..l :0. 00 (100.90,! 30.00091... .0»? 0.05.5103... 0.093.002... 0:0“.8005001 .0: . 3 A. ’0. V 04.10I‘ '0..-A .l A: . {...-1.0.0.0 . 0).... 02.00.30. ... . . . I... ... - .2005” .. ... 5%»...20 ...-...nufibuw fififiuwmflu. .Jsmfifiu.’ ...). .231... .3” Ls......_..a.....r....QA 3 ... :5... “4...... ...... ...... .. 0H0 ...«0..... 001.... 3.011 V... [0‘0100n40u ...h.u00}0900?33&u01 UV! . ......_......2 .2... .933 ... ...x......9$.r...c 1w... 0 . .... . 0 . . . . 0 00' 0. .. .0 I7 ‘3.... ....nnsma q.fi§.knlxou.%aon 0 .0 0. . 0 .0. . ...su. “33.00““ 20‘3““: haw L ”flak .00 " . . . ..0. I .0000... 0.0.0.? . . . . I . . . . . v0.V . .. .00: 000. .0-0. . . A. .0 .00 90.”: 0 0000.00,. .00 .0 9"..0 .10 .I0. . 0100...... 0. 3.. . 6. Ink Av 0 0 O ‘ 8 .I.a..0. 7' 0.340.013... .000 3.007... 1.0 0n..‘0.00.00tx0~‘.1. 00000.0 ’00.. . .0 .0, '3.‘0.‘..V>Alo .. W00 .0 9.00. 00 .000, .‘r 2% 1—03.04; 0‘30”" 00.00.. 0.30. . ,0... 0......‘f0... E000 0050.00.01”: . 300.... a. ..0.0JV‘.'.0|.009 00300000.. 00000 .. .0... .000». . 0‘0 .. ...... 0.... . .V.-00.8.5.0 .030... I. 0. 3...: 0 0.03.3. .. .2042. .. ....A 1.0.10.0)... 2: 12:03:33.0}. ‘ v 0.0.00'10 .007"; ’4’00'00030... 43020000 0.00 (.0. .000). ‘:3‘§30J0.L01.00002A00 o. 0.0“w 0006.00.00.10. . x. lust}?! {fit}! 5.3.30... .. "0... I! 0 o. . .....0 .... .33.!!! 0109.. 32.0.0031}: ....?:......... .......att............0..... 5.3:... ... .0 3...... . 35.30.0130,... n 0.00 03.0000... .0.0.000..007 .3. 0.0 00 004.2‘ 000. .1 ...‘ob.00..03$..~00.0 . ‘0 .‘0.’ ‘0 ‘. 0.0..‘r‘ 0.9,.00'00 0. 0.00 0.0.00 0‘0. 100‘. 0.000.0‘00 .00... V 0.... 0 00. . Z .. V .01... .. ...... .. u A A 0 6.3.0324... : 30.000.33.03... 0.0.3.0.... ..r ..V010C....€3 . 0...0..00ofi..00000.-.0..vr...§:.... 00... .l 0 . .. 0.4.13.0... .. 00.1.0.1... ......81 ......7; ... . .. . ...». .... ..0 . ...: ...... ......0. .....33 . A . . A .I 320...... . .00 "n...- .0C0000'O¢. .0033.-. 300.2020... 20 .. .1. . 0.0 $0.. 0 .. 0A.0.f...00.0:!..00.0ol.... 2..-..- 0 . ....00. $71... .A 0.. 3300...? .0 32:00:02.... ... . .... . . 230-1002111301. 0...... 0... V .. a. 00.003... . 0... .... 10.70... 0 .0100??? l...0 .....21..\..100 0.x. 9* ...-u. 000““! 5!...‘au0flo. . 0?... . .3. . 00. . .... . 0. . V. .. ...... 20.... ...: 000.1\.4‘.0.100 3.2.0.3. w. . .‘1...010.V...0 .30.. ... 220:... 30 . (3.000....zut... . A2... . .... .. . . V. . . 0...: 1.00 0.0“..030300 a, ..0 .E 30.3.2011 .0005; 3.030132152ull‘013. Alon. . ...... , . . .. A . . ...00,.._0..:.2:..Z. 4.0303001. Ah: 1.01.1.0 .... i=0. '0... V.3.a0«0000vr.0‘0b!0.!00..!0.0¢l..0.00 0 . 000.302.060.010. . . J: . . I 0. ... .V.:v’fi0$|\ I0.‘..3.I0 00 .000.0§0.0.01.J..0000.‘000..0.0. 0009.0 . «r I . . .V . ...Vét.‘ ‘IVQ. . “1.003.‘u.,. 30 .0 . 10000-00000 00....100 3w! ‘0. 3-0....\ 000... 90.... V . ......S.’ . 0‘0. . I...00..¢0.0..0l.0....0 0. :0...9.......00.§00‘0p0 000.0... 04.100000 0.30.0.1??? . ......0 0...: 0.. 50.93.003.030 000 00.1.0.0. 0 ‘2‘... ‘0? $3000.00 000.0 005.. 0.1%. ...: . .33.. .)1.n._40.\s"...00.10. . 0 1.700 . {0.20.01 3......1lscss‘ 0.... 0.0" ‘02900. ... I. 0.. A.“ d 60.23:} .V. 00.0. 0 0.00 l..\!.‘0u0. 3‘00...0.0V0I...a\........ 00’: 0.00. V .5000. ...: 0.: 0.00 ...0....$~0..& ..v . V 1,000.0. n ..0000.0‘.A.. O 0.0. .... 0.....0...t.2.‘..11.....00... V10 0 ..A\. . 0 07. . V .. 000000 finialh‘ ......0200 3000...}.340100000.‘ . ....1‘....fl.3 000.-..)000...003¢ . . 13.0.11..21..0M... .5. . ...-8.00.0.1 300020... . 0..0.01.00'..A00:00 300.00“ g.“ ...-00...... .2800: 22:00.... w... 0. 29...}: 00.0 .3“... A . 0 .. . . . . . . . . . . .. . . . . . .. . V .. . . . ...-‘0‘! .I00. 0 .0.‘ 770103.030. .....0...... ...“. . ‘0....l.l0...0;100 .0. .....00‘00‘0J0! n... .0. .3530ri000 .300 q 9. 0....00 . 0.0.0... 0.0.0.. 50010:, A .... .00...Y.0.40.. 0.0 .. 00.0 7.70 3.000.... 00.3 ..00....0.00.0..000...0IV. 0.0.0... 1.3.1:. . {00:0'5‘03 .000 05.,0‘0000. ...}... 000.00. ..00 0900 0. 00. 00.0300 0.0.0.0.. 000 .00. ‘01.... . 10.0 0 00... 0.00: ...0 A logic . I 0%‘09-d'0l0. 3:000: 00.0.0; 00‘“; ‘..0.0 0.10 0 9000...”.03. '0... 0.0.00.0.00 .-.0‘2000: 000.00 00' . A0... 30.0 o. ..000 .0 00. ..“»A$.‘.h 35.00.40”?! I...’00...‘Il0|l........0|..A.0.0-.I.0J 00..».0..0.‘3.. . 0 .0 .0. ..00001...\. , .. a. 03.0 .....I . .....0 . 0. . , . .. 00. V V .03 . . . A 0‘0... ..000000 0..."..03. 0. .. .UWu”0“0"l000"000‘fi.01“0.“-‘.0§"' 0 0‘0.0Qh00’ml«000 . . v 0 0 .000.- 0- 0V . . .0..0.. .. 09.00.”: .. ...‘0. A .00 .00.;0..000 .200; .00...‘.0.. ‘0‘}.0‘ .0... .0 100.003....70100 r0 0 $.00.._A..9....000.0.. . ...» ..0. . . .0 A ..f ..0. 0 0.. . . ...... ... 0-5.1.3003: 0.23.1.0 .V ......0... . ....l 00.03.000.00 0“ , . . . O.......Q¢0..Y.0...00 V3.10; 0A0. . .0 .0. .0. 00 3.0.000. .A ‘0. 00.000. . .. . 0:0 0 . 0.21.3.1. 0.00.0V01. .0 . 0.1 .t"? F. 0 ..0. .30.... 3300......‘00... I... . A \8: 0.... 1.2.0.9... 0. 5.0.2: ...... .... f... 0r .. . .0318...) . 2...: 0.005 . 10f}! .. 000.000.0000.. , .0‘0...‘ 3.000033. . . . . . . . . A. . . .....001 0:02 000.10: 0 .... . .V. .. $000.00. .A .0 . .1. .00.. .-.0. ...I. 2.0.0.... 0.000.C0 00. ‘00).‘005 0. A I“. ... 0:0. .900‘1.000300.800‘l.0'0 .0-003..V...[0. 40.00-09.50. . n . . . V . . V . .00.. .0' .0.!!.00(.............1.. ... .90... .0 \. 1 ... ... ...}... 9.... . hi}, . ' 0”. m I? 1’ 7072:3820... 2.1.0.0....00000.‘ '0 3.0.0....00... V . . .... . ...-000 01.0.0.20... .0... .00.. ... 0. .... .. 0" “0.0... 000 . ..0..’ .10. u.£‘..092‘10 1|). 0 0 0.‘.\$$‘ ’. 0.0 A . _ 0..q0.\«00'. 00000503000, 00 .0. 000.0.00010 . . . £0.00... 0.. 00.00.... ..090.u.00000.00._ 1.0.. 0.. 0 0 0 0. . 0000.. 0‘00. . “0. ‘00,0{ .00.. 0 V." .‘,‘ . I . . 1.0010. . 0.00 .... 0 00 0.000. ..0 7.00.1.2... ’00,.00. ,. u 0§0! 00.. . 3 ‘ '0‘” A. V . . 1.0.0 . ..v ...... 0.00. ._ .... 0.0.1:: .... 1. f . . . . . .0”, .- 3 i ‘0 0‘00 .001 00.0.04. .12....‘0300. .000... 02.0. .01....1090. 020.....0. 0.3.. . A . A . . . . .th-n"h0$0. .00 Juhwn..0...0.0,’0.'0’0.0 ‘ g A .. .0 .03..0.. ,. . .. . .. .. . . 0 . .[0 ...: .0... ...... . a. . . . .. . .. . . 030.0... .0001. k 0...»... V. . . . :00 30. 0 ‘2‘ 000000ha0ofi A , .5...» .V.. 20.00 4.. 000-....000... ..0.. 1 .010. .0000...’0.. 00.00.032.00. ...-.00 ... .’ 00 .70.}0’19 ..0’0 0 .‘C‘ 00 .0 . r u. 0.5!...1000 3.00800 .... :39... 0, 03090 0 0‘}. . €073.510IC .0 0.0... ... {03005.0 0 .0! A 0. .. 01‘3000‘01...0!.0\I:! 1'0’0030000‘0000'1‘0005 .1310. (.00... 00‘ .0 .11"... 0.0.0 0000 000.000.00.000. 0 9 0.00.00.03.00 0.00.0013? . . A . ... . . V. . 0.,00.0V.0‘.0.Q.u 0. . 0.0. 0-, .0: 000 I .010004.000Ir00" ... 0.0.0..10.0... 0"0‘00“ 0.0“}4 00"}...‘00 ... I.” 0... 01 .030... ...; 30 0003.00... .. . 030... 3.4000 00.0 4 A300003000...‘ J00.- {00 00......- 0.in A o . .000... 0§’-.0- 0.0. 0.. .0.2...0..0. ... 00.0.3.1: . . 02.00....Ol 0.0. .300 0.000.011.8030‘0300201100 ....A .00. 000000. V0000. .19.... 0. 000 00.0.0100 0. V0. .....000. 00:20 1310!..P5. 1.0008 . . ...} 3-.....00 0 0040.30.03.0300 000 300000;... .00" .0 ......30, .. 0001. .0 . 3'3; 1:50...- 1000090000010505‘.» 01.00 . . 0.300.973... .00. .. .0 C ...; 0.! 000200000000an00000 ixi’inax0ol’00 .0‘ .00110000. . . . 000.0 0 A A 2‘. 0 '60' .. 0. 1 x... 00F....0...0..:.f..00.0.....8I.1.00t...05.0.3.100333 . . . . . V . . . . ..l . 90... I... 39.25.150.00. 0.0. .. .. 4.010.220. . . A . . . . . . . V , ‘.§ .00. A . .. 02 I ......1...0..2v001.0«.30.: Spurn”... V . . . . . . . . . .. . . . V.A . .0: ..4 :05: ... V . ...Ihuuhuswtoflu..0.‘|.ut.i:.0.0..nuvl .04...0.0lo....0. .. . . . . . . . . . . .0 107.\000.0..0Y......00.... 30.3.0.1)... 212.003.100.02 3.0.0 030.31.“.QI..'10..‘ . vol\0.00030 .— 0 . . . .. .. 0‘40! 00 v.|0u‘... .0 )0 .0. .V.0~.0..0..0400’J. 0.000-000. .20AI’0 00.0.00“ "I. 0.0 00" .iigg A . A . . V.. .. 1.927... ....0... ...03........0.0.. 2? 33.35.30... .0]... a...k€0..\03... 0.... .- ..}...00. .00“ . . . A , . V . A .. . ........0... 3.0.... A900......._.1... 0 A. .....h... .0. . 193......0': 90.00.. 00......0. .e .....1 .31... 9..» . 0.0%..19....:? 2.0“...A....0.0.:0.01.:, 21.1.0.0...2..0.|. 9“.-.) 1.020.000.” 1.10.00.08.02201. 1...... V . A . ....0...7‘0...V . . 10.0 ‘20.... . "0,.0000~.‘~l00.0. .0 ‘0'}...r) 0 Oil‘01‘00 31). V0.0 .. 01.100. .0005 ‘30... 00......0. 0.00. 10:... $5L‘0A'.31000.‘00 1 00‘. 0... 1.000.903.000111. ...-0. .10.... . 3"} .30."...0. 1.000”. \d. . .00.! . '00 0. 0. .00 ‘Q 0,0003‘0 .00. 4 .0 .00 V 00 00 ‘0. .. 0. 0 .. . . 0 . . 00 0 ~0 00.. .‘0‘00’.1. Ito.1.47‘!0...0.00’\00..000. 3 nflmfiu. 00.00030001200000. , V : A . A . . ., . . 20.90.13.132? 0.0... . . . A .A V V 00.00. .. 0 .0. .01 .. 1 . . . 03000 0.0 0A..§.‘|...0000.. A... . , . V . . ...... ..0..,..0..002010......A. 30.0.00... . A 0.10.12.00.31 .. .0! 3.0000 .130. 7.00. . .90.]...TAin‘ 051“..."qv‘000020... .00£.0 v 4‘0 :03 00.000 .00!.§.0~0 . 1.000 .039... .. 0.30563 3.1.. .0... . E . a 0 .0001... .. 0 00.000 . 0” 5.00. ._ . V A. is i...?2§§3? 0.3. 0001 1.3.1.0010“ . 00.0000..AV.3900 ...0050..00J- 0 :IJCI.‘ :4 0000..v0....000§’00‘:00..00 ..!A‘04.10.'0040 0‘0.‘!\0<'000000"0£0104000 100001.000: ‘ 0 00 .0931» .00....) 7.00.. .30.... 6... 1.000 ..0..l.0........... 0 .0 (...-«3“...8.‘ . . A 02000000....0CQI‘. V A . . A A . . . . . 0 v 7.0.30.0}...0130. ... o ‘. i9ul¢ui.0!h««.0¢00~..0’00.nflw .. .0. ”10.0.0 4 1.0 040. . 0 ' 000! V . 0. . . 0 . . V 0. 50.0.. v ...: 00 . . .000 0u‘..,.0000000.§000 7"... ’0.§0000..0“0'0 w 0“. . A . ‘ .0‘000000000’ .. ..00... i 00005.. 3.0.10.0. Six-.500: ’00.:0‘0: r0600. 0 1 . I 0 0.00.0903. . 100.000.00‘00’00’3030010 l 04000 0“.0000.:l .0300’03V ‘0..V0011000.0.‘...':000dl‘.00 . owl-0000'. . \ ..la...v.§3300.000 0000.5...00000...00..00.0.0V.00.'00“. 31.0.5.0..013- 0.12.00.0‘0. .. . . V .0...1.V00..3.J..0...3|0.....00...0J.w..J....0r03.‘+... :0. 2.0.9... $1103.33: 0 0.... ... . . . 00.00.10.001... .04.]!000000000... 00“}..v0...0.. .. .... 0:03.......!..¥0000000000’00'..l.‘310.0. 0 .. y A . . . .. A [0000000000 0.0.0110. ..00 00.03.03.000 0.00 ‘10 .11.. 0 . . A . . . A A . .0. 0.700.720. .0 a. 20.00. . . -..... ... .. .v. ....$.20.........0.o..........V....«0A... 2‘s». 3.0000220}... {Swan}: .30 0.000 0 .. 0:00.. ... 0.0.0. .Q 00.9...0...\10.o.0..0§-..0 0.‘§§.0.0.10..... 00.9.0. 10’.- I. .. y 0. 2.;. ..-. 1.07.... 0 . . .. ...0...:........1.30....‘1‘021703: . 0. 0.0 0001030.}. ... .....0 )0; ..Vi..... . .424 3.1.00.0..01’43... 000 00..) $010.01“qu . V .00.!!0.....0......0.0..0..0 0.0.0000".- 2. . I 10.120.20.450 52...}... 000‘.0.01¢.0I.0..J§.~. ””000... ‘00. .60- ...’..Q05 300000.72... u.-‘0.00....0.0 ..009.‘ 01020tl'la 00.03000 ‘ A '00,. .\ 3.... J. 710.002..- .60. 3.0.24.0!!!)3. .000! . (‘30:. 0.0. ‘...0“00..0v.00.000§‘0.._0000000. 0 0‘ 1 0" 5‘ . 0.10.1 J,.,.0.8....900 .0. 4.0000. 00 .0000. 00.: .0 . . .0000‘.00\..!0000‘0\C0‘. g0. .....V..V . 0 I ....0 00"0,‘000. 0 ..0... 00 00.10033... ‘00). ~0’0.000.uv0o0.!..0.0..v .0: 0.. 30.0..-A00A000‘0.0.00000.0.0¢000-I.0000 000v .. 0.... . 00.0 0. 10000000 0000 A... 0 .. .0 .. 3.0.30...3...0.v).fl..uu.!..303. ... 0.... r -w. .ll000 30.13.00. .. 304:0 . 000... .200. '0‘ .0 00.000. 1000 000 31......0 ‘...0000.‘.00000'0-.00000...0000\I........3I. 00' .0. 000.3000 50?;0‘... 00000...‘00000I.000..0000 $53.01... 0:! .00....0...0_.... 0 .0..0.....0..01.1..§ 01...: ..h 02...)...010. ...? .. z 0:. ......V. .0V .......A . ... ...0....J0...1..\......0....0.; ... .. ... 1.3018... 00...}...0 .330. , . A 0.200.: ... "T “2.000.000: 0. 0.05:0’00_00.00Q V0300 ...4 0000.0 05.003.03.01...‘ . .3 .. ...! S. ...... taJ .30... .40.! 0033300....V..J:.0 $1.0)..9080 t... 0. 0.... .A V . .1... —. 001‘ 000.00.0‘ 00.. ‘0. . 001.. 00' 00 . 0‘00 0 .0000) 00.00 . 00s .0 00.0 00.0 .0 0. .. . v 7. 0‘..”.0.000 0.... 9.01000 0 ......0..000.lu.0.l ..0 20...)..3. 0.1.3200 infamy“? 0“”...2 0H3£0P . ......0. .....0: 1502‘.0...0..:000..\ 0.4.3.0....0001, ...! 00003.0..0000 0....000300010003’009-0300"- '04.: v0 3.003 5.0,......00...0l0000.0.0100.00.. 9.0.0‘003 ~...0.!...\0 000w. tillfiil’l‘u 440.0 "000$..«n00.0. 0‘0”... 0.0.0 .0 I .J..2l.000:00...0. .00.JJ.J.~ 1.0?) 1.0000. \Q..000I.100A.000000.0...80.}. 0. .0.....007.O0000.00‘00A. ...... (00.000000321205000 I0..004.:-00.0. 30000001000 0.0.0. . ..00.000.0.Q..\ £10.. 1.. 0.. . ...0000000000. 50’ 0.0.3 ‘ .. $.‘0.0030.0. .231»... .38)....Ill3; 3 ...-1...... . l).3.....J.0S0..-.9..000.1J.03.00 .0 ...-93:00... 0.0.0.0030... .... .0 ..V -C0.-.. 0 0. . 0.. ..0. .30.... 0.0 A .0 20..1........-..J.00i0.0..;.. 53.23%.1033 .... 0. "“1? ....00.......0...0..00...0..0§ 10.03... 33.33 ... ... 210...... . . . . . . . 0.0.» .......l0..010...0.......\ ‘0‘ 5.0 .J 3.4.2.1000. '00:: ‘0. 0.10:0 00000.0 21.3. , . A. . A . . .. . . . . V A .000...03~\.0.\..I0‘20010.0 030 000030000 .0.,‘,0oq‘o:vv0a"l 0.0 :0 ...-000. 0 .. A . .. A . .. . . €60.03. 2.00.0 ‘22, 0.003....100109033' 0. l‘§0-.0’1'0~..‘ ’1! . . . . . . . . . .V.-0.. .. ‘2 .v 30.1.0.0...»0301.30.01.53}: ...-.0 (0‘. 00.02.00.103 . . . .. . .0 0’ ..A: 000.\\00 ...v... ...10...’ .000.00i\0...0000‘0.1.r000.0¢‘.000. .04 0,000 10. ...}..0‘000. . ...-0.: . 000.303... ...-‘00.}...010 0hr.-. .000'. 0.00.10.20.01 ....10.~.i030..l3r0}00 0.003.300.0013... 0 :2. ....\ .\\.l.0....0..)..\1(0 I .00 "00...’ 00. 0.0.0. l.0.95‘53‘.0.1.n' .i‘.0“.0h00“00u0|0. .00 . 0.300.190.332 .1? 0......200.0...$ 00.. .....0‘3..’0.00$.Jla .00 . V 02.00.01.103? ......0. 4 3070.0133053002. «00.0..0‘2..00000§00I000004.0 .0)L % I: 3010...; ..A....11.:.0.... ..A 0.......S.0.$ .... 00.1.“ .301......l.$..l.00.00.000 .. . . ... 1:... ...!Sqa...\..0- .7100..\.. .th...:a........€0.0301 .00.-"‘01 $0000. '1‘00‘0000 00100.00. ‘00. ......0.‘ 30:00.? ...-.30. 0000......Y0100 000‘. 0) .0 . .0. .0 . .l‘. . . . _ . J. ....1‘0..00........21)0.2.2.203393100..0. 2100.... .000 30000.. 0.0.0:... . 0.1.x... ...I..!3.0 1.3.0.3.... 3.0.0.! 03.30.0002..£.,.0..3 ....0........:.0.00. ..I....30...3.0 ..l... 0 0 1.07 01.000.20.23... ... 0...‘00.A0... .. 1.0.0.000 00".V...0..1.. 000.0'70 200.00.001.‘ 0.30000 vii-V 0‘000 O..0.00~.l.’3\0000.. . 00...: 0 0. A‘ 0. . .010 0.0 .l‘! A. .0 ... 10.00.00.005 . r 100..., #:0000000: ......00 0 0-00 00-001000000010.00..’ 23‘ 90.33.3043 '00 V . . . 0 301040.. 0000...]...- ...00! ..0’.—.l....r 0070.0 0000...‘...09.0’0000. 700’ .0. .0. 00.000. 0.0 a. . . . . . 0|. '0 0.1‘000)0000.000= ..0‘00)!.... 0 . . 00095.“: . . 3.02.0”. ‘0... 0 0'0 .0 .V . 0.50000 4.0 0. . 0.32.0010.” V0...00.l..0..s0.‘.v_ 0-0. .. 00 00!... .. 0000 0.00. A ... .00 .0000..- ..0.... . I v ..00.....t..\...0. ...... ..0. . . .0 , 0. . . .. . . ... . 0 . Q. 0000 . 0.00150. 0.0 0.0.0.0.00.’. 00. 00.1.0 0...! 0....0‘3..-0.0.-0.0. 00.00. is 0 .0.A.010.0...0.‘0 1.. ..0. . 0‘ A. .000! 900‘“0:...'.4’..‘0.Q.‘00..3 ,‘Wg0; \0\.\. . ...? .0...00......A.0. 0.0.0.. . ....0......0?0A.A. 000 5.0.0.0.... ..0.05..0 01.5.7.0.- .0’0..h.§‘0030.000. ...10‘1100‘.:N.0§0002300 . 0. . . . . 1,0. 19,000.. . .... A. .0... . 0.0.9.0 . .00.. 0.0%, 000.. 00": 09.00 2‘.‘.a‘!§0" ‘03‘égg 00 d..00.$......0. ..000. ..., 0 . 000..., 0.1...0577‘0‘500000 .2! 30...... 3.0.0 3’. 2 . . . .0015 0:00.... . 0.... I....0...0~i\1.‘.00 1/0 :90 ‘0. ‘i..09.v0...0 i‘....0. 1.10.5213... iii-3.1.50..0A. .70....5 300.01.300.08...- ‘10...~ . 00100.0 .000..0...A.J00.000100.0~00 0.0:". 00.1.0530 .. . .05.). .e 7.00 0, 0.0’3’0 0:000“ . .00030000. 10000. .0. 0 0 00 00.. ..0 0’. 000...:1 {.0 0.000 .0.00‘:0.O0.0. A 00.0.. .3030 ‘8.0....u.lt. 1.0..9 03.012.70130...‘ ... 31.0.0.1 0.. V . . . 00.0. 0000000 .0030 .3400. ’00330 0 0 0 000- . 000.: 0.0.0»: . 0 ~ ..00 $0000. 00 00.. 0 0.0 A 905.0..00‘000.‘ .... .000... 000.0 .. 0... 0.3.00... . ... .07 . A, . . ...... .mtu30n70233 000 .500. 0000‘ 0000.000 0. ’00-..00‘0 '00.. ..v|.0000100000 .... 000. 0‘. 0 $0.0. 5000.00.0000.’000 0 00 (:00. 0.0. .0}.I 00.00 000.0 0 000”00..‘0 ’i‘03300033 00-...000Ai ‘00..00.00..0-..O 9’ .000 0 .000.- 0._..-.0.I.0ll.. 0.0 0) .0..0I-r0000.-00..-, .100. 30000310...) 0......00...‘ 00.004.003.0‘.00000000.’ 0000.03 03000004000006.» A . ,0 . 10“.... 103.200... .0000 .01.. 2:030}. 0000.. 0.30 0 n 0 1 ..I .01’0090...l...0’.0..6000.0.00 00.5.)Xi00uzu0 0.. 000‘.‘ ...... .00... ...-00,010. .00.. 000.0...000.0.0.0 0 I 0 ... 200.0 0.0000... .’0. O 0.0}A...\.’.\Q0.0’0. .40 0.. 00.000030010. I 00:00. ’03:..00000‘O. 0.... 000.. I. {...-000.00.00.00 ‘30....00 .. 0 90.10.100.0000 0.00 0000.000 .00.... .000...) - I..0s0o000.0l...09:0.0 0000.03. 00... v.00. .... 570.000 . 0 .. litre 0.0.}. 00’000000.‘0.0¢.00.00\0‘....““0100 . 000.300.003.00200000 0 01.1...‘000. 0:10.00 00.0. 30.000 . 1:0: 00......0....0_..00 ...-0.10040. A . 0000000003.... 0002.00. ..00 ’00:}... ,0 00 .0. . A Q. .03 .50 ..0030.‘..0.....0000!.J.0. g0...‘.0.4..J.-00...a.000:?0.. J00- ..Q.".Jnuw.h0.01~.. . , . . 30000005.! .0“ 9.0.‘00031l . V . . V . . 0-00... 0. 0 gt}: '5’... .. . .. . I... 000 000.30.000. . 0. 0 000A ....0 .V.... 0000 0 20,0 ‘0 ‘00:: .2 00.0.. 5.0. 0 0 00...}. 9. 0 A . ..0. QCPAC “.013; 9000-! 0’00 .. .. .0 ..0 ... 02.. 023.0: 0. 700.00..“ 0h, .. ..." 0. 0 . 00.0... . .A. . ...-......fi. 0 . 90. .0. I 00.. ...!0 040.. to... v _ . . ALI . v. 1.0 3 0.010. . V 000.30.11.03 0"“..0 0.51. ......I ...... ... z . . A . . . . 33.20.00“... 4&0?! . .. .. . gill: €00‘0.‘0130.s . . . .. . A 0.0.5.0...10.0000000....70 ‘03)] 010.20; 0-. .01....0.‘..0.‘.003..000§0 . 0. 0.0.0.0003 1000.500 .0‘.‘ ...: 0.200020..‘... 0. r... 00 .00 .. 00... r033 0... .0 A00 .0 m l 00. O .0 '00 00’IA0, 0 I ‘9. 0.00.".0..‘0. d. A. .. . . . ‘.0. I 0 .0 0.000 0..0000.0000.0 . 0. 000 00. . . A . .. 39.002.301 . . . A V . .. 00. .0 .00 0....- ~ 0 a 2 0t . A. {00. ‘ 00 .. . . ...; . A. . .0... a ..su .u0’ L .unn’.\‘0...0.1‘|00.-II \10 r...002§93.‘-.|.§Q¢ 0...... 000.931.... .0 .01.... .'0D 0‘.,\ U 000 . 0.0.0.0... ... 0. 000.000.0000.... 000.. A .00 03.0.0.0 ... 00. lo ..000 0 .0 00000012000000. 00- V’Q70.0.00..u000.00.00 0 00...... 00.5000. 0.. 00. ‘.0003.000.0O.000000....00. 1.0.0.30 .0,!€-0!.I 00". :080. v '0 ..0000‘00.0 . ‘3; .0... A < 00 01.0....- ‘0. .... . .QO‘O..‘3.9. A. o: 0...! . . 00‘. .. ’33....v. 0 ...00000.: 0 p . . DI I. 0 . OA 0000300.}: .0 '00 .0. 000 . ..00. 0000.00.00. . 00000.0(... 0 s .00.. ... .‘0000.00O\. .00 A ..00000. ..0 -0... 07.0:"000000’ .00 §£00:0t00 .A 00.0. 0.. ..n..:,00....\ 00.500005... . 0 .0000 00.30.. 0 I I0.’ 3:...“ . 000.00.A0....0o. 00. ’0 0 0.0‘\ 000-0000l.“ . 0:. . V ..nhw.a..$...000.£....a .s.\fi0u.fl.m:.:.‘II| ~ ...004..0.0..0...§‘ 300.01.]. 121310}- 10.0...)02...0...0V.V0.0‘.‘001..t . ....01 . .00...000 .‘ I 32.1.0700... 0.0- ......‘k 0. ..0 .. . . ...000 ...:Iag .0. 0000.10 00.0.30... . V . . . . .. . .. .. . . . . 0.. 320.1003. 000 001 I A... 3.3000000? . A . . .. . . . 00..., 01070. . 00.0 .00 V 0.0.000. , . . .. A. . . . . . 1.000;)...03503.‘ ......20, ...-.... 0 . . . . 0.03.1.0}... .3000}! 0.0.7.000 0.0.. ‘0. 000...... .0 0.0. :01; ..‘00..0 0-0.... 0 0.. .00. 60.5..- .100. A .0. 0 0.0 0.101.032.07300 . ... 0 . .. 0 0 ‘.00.0‘O.I; '00.. .00t00'\ 000700.000. . .000 .0 0.-.!013 0300‘ .. 0‘ .0010 ’x)‘. .2000... 0300.0. 00.0... .0...- 0 0.. 0. .0 .0.. 0000s.:0 .0 I’. .0. V’OO 3000.00 .0 000‘. . .3000 01000‘00.. '00 .9‘00: 0000 00 0. 1.00 000 01... 0 1 7:0 . 0.‘.. 00 0.1:.‘3. to 0.0. .- 50,0... . . ..v ... 0000 020.00."! A u........0.¢'\0.000 0. 0 .0 '93 .1. 00". .‘ K..000.0.00'000..0.§......00 00.05.00.,.000...000.0i.0..0...0..0 0 0 I .0000... 003 .00.. .0. .9000.....'.. \0 00 .‘o...r. .OVJV.‘0. .00 to, 0.0.0:...00 0 \0 .... . 00.0.00. ’.r00)... 0 . 0... .....‘.00.0’ 00000 . .0 . ... ...-0.. . ..0.0_0.....0-0."00a0.0 . .0 .10 0.0. 3.00.00.5000.’ 00.010...‘ ILA. 00.0 . A .. . .. A V .0000. V0. 0.0 . .. . .0... i... .V.. 0. 60.0.0-- A0..0.0...00.'.l.0; ...0. . . . . .. . 0 .00.: 0.0...3103. O u .0. ...v... ..0 .00...0.. -000 7.0. . . . V .0 ... .. 00.0.... . . . . . . 0.0....1. . . . I 0.0 0.‘ . 0 9'0: 0 ..,A. I. 00..400 O. 0 00.0, .-..0 0 0.0000. .....1’10 0 0.0... .. 00.0 ... 0. 300000000: 0 0 0 .00.. ...V00 0 . 0.‘0.0 .. . 0 00. ‘000501. ‘0- . V ,. ...50..1.A~..0...u .. . ....- 0 0110.00. 0.0 3.0.1.0019 35‘2'7 20 ..A.. 00 . ..s .... ..0 00.‘.0 0.0 0. 303.000.1000 ... .‘0‘010 ...A...... w 0$0....5..Iv!.0d:..;.‘0...\00 0 1 , . V. . . V. ... 0 00.00920. .0 00.. 0 .301: . 0.. 00. 0.0.0,. 0. . ..0.0 0 . 000 0... I. : 000.00 ‘0. .0000. 0090. .0; 0' ... . 0 o 000 0.0.0.0 .. .0070 . J 01... A. 0.. 0 . ~00. . . 0 $000 . ....0. o 0 v... 0 00. 0.0 .0.0.0 0. 0.0. . ..010... . 00.. V. 0.0.0... . 0000.. .0.. (.0 620000 000 0 0.0.0 0-. 000 0.0.0 A0 0.0.0.... ...A.:. In. .00 ....I 1 3.00.1.0AI0...‘0...000.‘0 00 .990 ‘0... .. . 0:00....- ...0..00....0.000 .. 006900....100...‘ . . A. 0.... V... .00... ... 0. ..‘..‘0\0u . 0..A ...0 0. ..0 ....l ....0. . 0 . . 0.0 .. .0 . 0.. A . 0. ... C 2.0.50. 00.301.012.110. 0.. I ... ...-0. '00 0.00.00... . . .. 0 on... .....0... 01 1.0.... ... A...\..A. .. 0. . 00.1.. . .. 00....- I . 08.00. . o . T... . .....I. 0 ... .. 00 .. ... .00. 0.3.0000 ..DJ. .H y .000 . . ... .A 0.2.... O0 0.0...‘00..000.0..0 ... ....0. .... .2.- 0 n . .0 0 .3050 a . 0. . A .....I- . . .V .0 .A..v0|.0.0.n..l0 .. .. 00.030000. 0.01.050). ‘0 . . A .A .000 . . . '0 0000.00 «07 .l. ‘6. 20.0! .0000 J. ... ...-0 L". .A .0. 00'. .I.'l.. . ...0 0000003.. 030! . 0 . 00 . A . .0 whtw....0n.‘“.0. .. I). ‘3. 0.00 . A . l....00 0 .0... ... ...: 0 .. . ..00 0 ... 0.0. .030. .0.r. 0. .0JA . . ... . 3.0..‘0. 0 ... 00. ..01. n 10.... p . .00. 0000].0\-v.00fi0.. .. .. .- 0... . . 0. . . ......0 n 0'00.‘00.... In, .0...00\. u .‘0 ... 0:! 0.000.? .0 . . . ... V A . .. ... .. .. .0 00.0... 0. .0. .000. .00. 000.- . . ..l . 50.. . . 03’.‘ .0 0 000 .V.. 0.000 . n0... 0. 0 .00.... .. I... .00 . .10... A . . .000. I . 0 . A ..100..0..00.000.‘.. . .. 0.0.0 .000 I... . . .. ..o... 0 10.... . . . . . . ...0...0.0I00. ... 0.0!...Fos0. ... 008.00). ..0‘. 5'0. ..00 . ......0. ...0 . . . . . A 0v..0.0.0| . ...,o... 00 ......00. 00. 0 C A00... .0 . 03,0109 .. .. 0.09 -0 ... ,0. 0 .. .00 l000-Q...O..00A0 ..0. . 0 0000.0 )0. ..000 0 0-0 .2 V, . 0.0 010 . 0.. 0.0.: V0. .. 0.0.00.3; ... 0C. . .V. .0... ..0000.-.0. 3.000000; 0.00. . . 00' ... . . . '0 0. . . .. . 00. . 0 .V v.. . .0 0 90.01” WI, . I O . . 0 .0l 0 5000.00.00. 0 .. 0. . .0 . 0.0021032. .0 ... 00: . .0... .....0. 0 . ... . . ... ...-.0. ..0 A 0.0.00..0‘IO\I..0. .100 .v I. .... . 0 ....I 05‘ '00 0.0 0. .. . . ..30.‘ ......00. 0.0 .0. V v 0. . . . :40... 0 .0. 0.0 000.’..0 . J. . . A . 0.0.3.0 .0 0 .10 0». ..V... 000 .20.». .. 0 . 0.... \l... 0... ...... .0 0030 ‘0 0. . . ..A . . .v . .00.... 0 .A .0 . .... 0.. . . ...0... ~ 1 0.0 .. ..000 . 0 0 . . .\. :0! 30.93.000.01... .0. . .. .I 0 0. o . . .- 3.2.00.1. .v 2. v0. 0 10.0.. . .0..0...0.0.0 0 .... .00. 0.... :07; 0 .... v _ 0.... ..0 0 ...A. 9’... 70000000033... 0 .0. .0. . .0. ..0 . ...'-.0.0..u0...5I00..L.0 0. 0 0030.0002900: 1 ‘).0‘30.0 0.30 0 0 . 0 . o u ‘0 r . ...0 0t . . . . . ... .. A . A ... 0.00. .. ... . I .‘o I ...: .... . . 0 0 040.. 0 3000.0...00... . . ...90 0. 0’00. ,0..I\§0....0... .. ...0.0.'09010 . ‘00 .0 . 0 ‘10. I . . . A . 0 .. . . . .....0..xno.00...0 .0~. . 0 \ ‘0 10.0.. . 0 O ..A..... A. :01‘ I..0I 30.0. .00 V ...00 . . A 0.700 . I . . .. 0 I . . . . , ..0 .V‘. 000.00.|.'0.JA. ..033. 0300‘ 030.001.039A ... 0:003.‘ . .0.0.A00...~....00 I0 ... I0........ .00.... I .' .. . .7-.. . . .... 0 . .v .. . . ... 0! n . l. . . ..0.u, 0..0 . .......l~ .... .‘ .. .A . 0. . . . . 00 . 00.0 .00 .....A ... 00. ..f 4 AV 10. 0 . ...-.r0. 03“.’.u.0..00...0 ..0 ...; ...NNU 0.100." “0.0.0.3061"? .. ...-00.. 3.....0...0.0.000.. "0.. '00 .0. ,Zlnol... .0 39.1.0 2. at... .. .20.... .. ... . 0...... . 0. .20.. .,. ...: {$01.30. . 0' . 0.. 00 . .\ ...... 0. - A ...-.. .0 . 0. .... . .‘t....0...0 0-.)0 0.3» 4......) ..mt... ...... .moltidfifiiflm 'J‘vwfl 0.0 . . 0 A. 200.... I n. I... 0!. . ......V , ll. .97.. ....0....0..;.. .00. :00- 02 0020.: . ... .- ....0...0.... ..0.... 0.0 I... . ... 0. 0.. ....f. . 220;. . ..0 .... ‘0.1 .A 0 . 0.0.50 . .00.. 31 I '0 0.30.0.0 0290.00.00.01... 0, 0‘....0‘.0..\ £000. .0.. V. .124». .0“ 0‘. Du. I -0. a . u 0“ 34‘ . . . . .. . .00 0 . 0.. ......000 . ., . _ A 0 . .. . 000 . . 0 . a .00.- . . 0 0 0 .. . ...0. .. . . 0 .1. . ... .0 . Q... . I V 0. :0 .0 1.... 0 0 9 O... V. 00.... 0.00. \ 00.1.00 . 0.0.. Q; ‘0. u . 0 ‘00 0 0,0 In 0 . o 7 . 0 . . ....10 0 0 I. .0. 0. . 00.. .. .. . .. 0 0 I. .0 .90 0. ... .. 0.0.. V0 .000I...0o....1 .'..A0..000. ... . A . u . . . 0 . .0 o0_00..;\ . _ . . s. .00“.. "0‘0 .01- ....... .... LIBRARY Michigan State University This is to certify that the dissertation entitled -.-c-o-I-I- CHARACTERIZATION AND PROGRESSIVE DAMAGE ANALYSIS OF QUASl—THREE-DIMENSIONAL COMPOSITES presented by LIANGKAI MA has been accepted towards fulfillment of the requirements for the Ph.D. degree in Mechanical Engineering M of; Major Professor’s Signature 7/25// 0 Date MSU is an Affirmative Action/Equal Opportunity Employer PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DAIEDUE DAIEDUE DATEDUE 5/08 KIProj/Acc8Pms/CIRC/DateDueJndd ——_— CHARACTERIZATION AND PROGRESSIVE DAMAGE ANALYSIS OF QUASI-THREE-DIMENSION AL COMPOSITES By Liangkai Ma A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Mechanical Engineering 2010 ABSTRACT CHARACTERIZATION AND PROGRESSIVE DAMAGE ANALYSIS OF QUASI- THREE-DIMENS IONAL COMPOSITES By Liangkai Ma The applications of fiber-reinforced polymer composites are gaining accelerated growth in aerospace, automotive, and emerging energy industries. This dissertation discusses issues concerning the characterization and progressive damage analysis of an innovative quasi-three-dimensional (Q3D) composite using finite element methods. Firstly, the average field homogenization method based on finite element analysis of representative unit cells is implemented in ABAQUS/Standard to compute the effective stiffness of the Q3D woven composite and compare with the conventional laminated and woven counterparts. Particular attention is given to the optimal choice of geometric parameters for achieving high in-plane stiffness in the Q3D woven designs. Secondly, double cantilever beam (DCB) and end-notch flexure (ENF) tests are simulated using the cohesive zone model to investigate the progressive inter-laminar delamination growth through the undulation regions in woven composites. The effects of undulation on delamination propagation in woven designs are highlighted. Next, a three-dimensional continuum damage mechanics model for the prediction of the initiation and evolution of intra-laminar damage mechanisms is developed. Combining this damage mechanics model with the cohesive zone model, different damage modes such as delamination patterns and fiber-bridging in laminated composites subjected to transverse static and low-velocity impact loadings are simulated and compared with experimental results. The necessity of in-ply matrix cracking is emphasized for the accurate prediction of delamination in laminated plates subjected to transverse static loading. Finally, using the developed three-dimensional continuum damage mechanics model and the shell-solid modeling technique, laminated, two-dimensional woven, and Q3D woven composites subjected to transverse low-velocity impact are simulated in ABAQUS/Explicit. The Q3D woven composite is shown to have higher impact damage resistance than its laminated and two-dimensional woven counterparts. DEDICATION This work is dedicated to my father, Chunming Ma (1943-2003). My father’s integrity, responsibility, and courage left an indelible impression on my life. I will be eternally grateful for his example. To my wife, Yiping Zheng, I thank you for your sacrificial and unconditional love. Without your understanding and support, this work is impossible. To my son, Guanghan Ma, my precious treasure. You bring me more joy every day. . Being your dad is the best reward I To my mother and parent-in-laws, I thank you very much for your sacrificial and unconditional support of all my endeavors. iv ACKNOWLEDGNIEN TS I would like to thank my advisor Dr. Dahsin Liu for his financial support and guidance of my research. Under your guidance, I have become a better researcher and presenter. I would also like to thank Dr. Ronald Averill for sharing with me his knowledge in the modeling and numerical analysis of composites. 1 am very grateful to Dr. J ung-Wuk Hong for sharing with me his knowledge in damage mechanics. I am also very thankful to Dr. Alfred Loos and Dr. Zhengfang Zhou for being on my committee to guide my research. I would like to thank Dr. Craig Gunn for reviewing this dissertation and making valuable suggestions. I owe thanks to Mr. Fred Hall and Dr. Dirk Joel Luchini Colbry for their support that makes the computational work possible. I am also very grateful to all other faculty, staff, and friends here at Michigan State University who have taught me new knowledge and helped me during my stay at MSU. I would like to thank Mr. Michael Dewald at Dassault Systémes for allowing me to use ABAQUS V6.9-EF1 and for his great support of the last part of my computational work so that I could finish the research on time. TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ xii LIST OF FIGURES - viii Chapter 1 Introduction - - - - - - - -- - -- - - 1 1.1 Background ......................................................................................................... 1 1.2 Statement of problem .......................................................................................... 2 1.3 Research approach and organization .................................................................. 3 REFERENCES ............................................................................................................... 6 Chapter 2 Characterization of quasi-three-dimensional woven fabric material.. 7 2.1 Abstract ............................................................................................................... 7 2.2 Introduction ......................................................................................................... 7 2.3 Homogenization literature Review ..................................................................... 8 2.4 The average field homogenization method ....................................................... 10 2.4.1 The repetitive unit cell .............................................................................. 10 2.4.2 The average strain theorem for periodic boundary conditions ................. 11 2.5 Implementation in ABAQUS/Standard ............................................................ 17 2.6 Model calibration using a plain weave problem ............................................... 19 2.7 More validation with experiments .................................................................... 21 2.7.1 The specimens ........................................................................................... 21 2.7.2 Material properties .................................................................................... 22 2.7.3 Results and discussions ............................................................................. 23' 2.8 Parametric study ................................................................................................ 24 2.8.1 The unit cells and material properties used ............................................... 24 2.8.2 Parametric study by varying a/ h ............................................................. 26 REFERENCES ............................................................................................................. 28 Chapter 3 Two-dimensional progressive delamination analysis of plain woven composites -- - - - - - -- -- - - -- - -- - - - 32 3.1 Abstract ............................................................................................................. 32 3.2 Introduction ....................................................................................................... 32 3.3 Literature review ............................................................................................... 33 3.3.1 Literature review on delamination ............................................................ 33 3.3.1.1 Fiber orientation effects on delamination resistance in laminate composites ............................................................................................................. 35 3.3.1.2 Fiber orientation effects on delamination resistance in woven composites ............................................................................................................. 35 3.3.2 Literature review on progressive damage analysis ................................... 36 3.4 Delamination modeling ..................................................................................... 38 3.4.1 Damage evolution law .............................................................................. 39 3.4.2 Material properties for the cohesive interface .......................................... 41 3.4.3 Delamination initiation and propagation criteria ...................................... 42 vi 3.4.3.1 Delamination initiation criterion ........................................................... 42 3.4.3.2 Delamination propagation criteria ........................................................ 43 3.5 Numerical issues of delamination analysis ....................................................... 45 3.6 Model calibration .............................................................................................. 46 3.6.1 The specimens and material properties ..................................................... 47 3.6.2 Selection of elements and mesh size prediction ....................................... 47 3.6.3 Comparison of results ............................................................................... 50 3.7 Selection of mixed-mode delamination propagation criterion .......................... 51 3.8 Case study ......................................................................................................... 53 3.8.1 DCB test of two-ply laminates .................................................................. 55 3.8.2 DCB test of four-ply laminates ................................................................. 57 3.8.3 Comparison of plain weaves with 4-ply laminates in DCB ...................... 57 3.8.4 Undulation effects on delamination propagation in plain woven ............. 59 3.8.4.1 DCB test results and discussions .......................................................... 60 3.8.4.2 ENF test results and discussions ........................................................... 61 3.8.4.3 Energy dissipation due to viscous regularization .................................. 63 3.8.4.4 ENF test results and discussions ........................................................... 66 3.9 Conclusions of two-dimensional delamination study ....................................... 68 REFERENCES ............................................................................................................. 70 Chapter 4 Three-dimensional static progressive damage study of laminates 74 4.1 Abstract ............................................................................................................. 74 4.2 Literature review and the modeling methodology ............................................ 74 4.2.1 Continuum damage mechanics models for composites ............................ 74 4.2.2 Progressive damage analysis of laminated plates subjected to transverse loading ...................................................................................................... 77 4.3 Three-dimensional progressive damage modeling methodology in this work. 78 4.4 Development of the continuum damage mechanics model for intra-laminar damage .......................................................................................................................... 79 4.4.1 Formulation of one-dimensional continuum damage model .................... 80 4.4.1.1 The elastic softening model .................................................................. 80 4.4.1.2 Mesh dependency and crack band model ............................................. 81 4.4.1.3 Damage variable calculation ................................................................. 83 4.4.2 Formulation of three-dimensional continuum damage model .................. 85 4.4.2.1 Damage modes in fiber-reinforced polymer composites ...................... 85 4.4.2.2 Constitutive model ................................................................................ 85 4.4.2.3 Energy dissipation ................................................................................. 88 4.4.2.4 Loading functions (damage initiation criteria) ..................................... 90 4.4.2.5 Damage activation functions and consistency conditions .................... 92 4.4.2.6 Damage evolution and damage variable computation .......................... 94 4.4.2.7 Viscous regularization of the damage variables ................................... 95 4.4.2.8 Solution algorithm and material tangent constitutive tensor ................ 96 4.4.2.9 Implementation of the material model .................................................. 97 4.4.2.10 Modeling of delamination ................................................................. 97 4.5 Calibration of computational model ................................................................. 99 4.6 Case study-impact induced damage modes .................................................... 103 4.6.1 Model I .................................................................................................... 103 vii 4.6.1.1 Results and discussions ....................................................................... 105 4.6.2 Model 11 ................................................................................................... 107 4.6.2.1 Results and discussions ....................................................................... 107 4.7 Stacking sequence effect on the delamination area ......................................... 110 4.8 Conclusions of three-dimensional damage analysis of laminates ................... 111 REFERENCES ............................................................................................................ 1 13 Chapter 5 Three-dimensional progressive damage analysis of woven composites subjected to impact loading using shell-solid coupling technique 115 5.1 Abstract ............................................................................................................ 1 15 5.2 Introduction ..................................................................................................... 1 15 5.3 Proposed Methodology .................................................................................... 1 18 5.3.1 Shell-solid coupling technique ................................................................ 1 18 5.3.2 Material modeling .................................................................................... 1 19 5.3.3 Solution technique ................................................................................... 120 5.3.4 Cohesive layer modeling ......................................................................... 120 5.3.5 Implementation of the material model ..................................................... 120 5.3.6 Modeling of delamination ....................................................................... 121 5.4 Model calibration ............................................................................................. 121 5.4.1 Material model calibration ....................................................................... 121 5.4.2 Dynamic modeling calibration ................................................................ 124 5.4.3 Shell-solid coupling model calibration........ ............................................ 127 5.5 Case study of impact- induced damage ........................... . ............................... 134 5.5.1 Comparison of laminate, plain weave, and Q3D weave specimens ........ 134 5.5.2 Result discussion ..................................................................................... 141 5.6 Discussion of the delamination initiation criteria ............................................ 146 5.7 Final comments of three-dimensional dynamic damage analysis ................... 149 Chapter 6 Conclusions and recommendations 153 6.1 Conclusions ..................................................................................................... 153 6.2 Recommendations ........................................................................................... 155 REFERENCES ............................................................................................................ 156 viii LIST OF TABLES Table 2.1 Plain weave unit cell geometry [22] ................................................................. 20 Table 2.2 Tow Material Properties [22] ............................................................................ 20 Table 2.3 Matrix Material Properties [22] ........................................................................ 20 Table 2.4 Comparison of homogenized material properties ............................................. 20 Table 2.5 Unit cell geometry [47] ..................................................................................... 22 Table 2.6 Tow material properties CYCOM 1003/W-490 E-glass .................................. 23 Table 2.7 Matrix material properties from Cytec Industries inc ....................................... 23 Table 2.8 Tow material properties used in the analysis (Interpolated) ............................. 23 Table 2.9 Comparison of results with experiments E11 =E 22 (GPa) ........................... 24 Table 2.10 Tow Material Properties [43] .......................................................................... 25 Table 2.11 Matrix Material Properties [43] ...................................................................... 25 Table 2.12 [O/90]s laminate effective moduli [43] ........................................................... 25 Table 2.13 Unit cell dimensions for varying a/ h ............................................................ 26 Table 3.1 Geometry of the specimens [31] ....................................................................... 47 Table 3.2 Material properties [31] .................................................................................... 47 Table 3.3 Summary of result comparison of model calibration ........................................ 51 Table 3.4 Summary of mixed-mode effect study in DCB of [0/90] laminate .................. 53 Table 3.5 Comparison of DCB results of 2 ply laminates ................................................ 56 Table 4.1 Material properties associated with each damage mode in three-dimensional model ......................................................................................................................... 94 Table 4.2 Material properties used for three-dimensional progressive damage analysis 100 Table 5.1 Homogenized material properties for global shell section of [0/90]s ............. 129 ix Table 5.2 Material properties used for local solid and cohesive sections ....................... 130 Table 5.3 Comparison of forces obtained at delamination length of 40mm ................... 130 Table 5.4 Homogenized material properties using tow material Table 2.6 ..................... 137 Table 5.5 Homogenized material properties using tow material Table 2.8 ..................... 137 LIST OF FIGURES Figure 1.1 Fiber-reinforced polymer-matrix composites studied in this work ................... 2 Figure 2.1 Composite structures and their unit cells ......................................................... 11 Figure 2.2 Implementation procedure of the homogenization method ............................. 18 Figure 2.3 Unit cell of plain weave ................................................................................... 19 Figure 2.4 Unit cells analyzed with 12 plies ..................................................................... 21 Figure 2.5 Specimens analyzed to match the dimensions in experiments ........................ 22 Figure 2.6 Unit cells used for parametric study ................................................................ 25 Figure 2.7 Unit cell for five harness woven ...................................................................... 25 Figure 2.8 Normalized in-plane effective Moduli by varying a/h .................................... 27 Figure 3.1 The fracture modes and test methods .............................................................. 34 Figure 3.2 Bilinear traction-separation laws ..................................................................... 40 Figure 3.3 Flowchart of progressive delamination analysis in ABAQUS/Standard ........ 43 Figure 3.4 Cohesive element mesh converge study using a DCB test .............................. 49 Figure 3.5 Stress distribution at the delamination front in a DCB test ............................. 50 Figure 3.6 Comparison of predicted and experimental results ......................................... 52 Figure 3.7 Models used for case study .............................................................................. 54 Figure 3.8 DCB load-displacement curves of 2 ply laminates ......................................... 56 Figure 3.9 Post-test shapes of 2 ply laminates in DCB ..................................................... 56 Figure 3.10 DCB load-displacement curves of laminates ................................................ 58 Figure 3.11 4-p1y plain weave models .............................................................................. 58 Figure 3.12 Comparison of DCB load-displacement curves of plain weaves and laminates ................................................................................................................................... 59 Figure 3.13 Three models of plain weave [90-O/O—90]s .................................................... 60 xi Figure 3.14 Comparison of DCB Load-displacement curves using the three models (1 and 2 denote the undulation locations in Figure 3.13 (a)) ............................................... 62 Figure 3.15 Comparison of ENF Load—displacement curves using the three models (“1” denotes the undulation location in Figure 3.13 (a)) .................................................. 62 Figure 3.16 Energy dissipation in DCB tests .................................................................... 64 Figure 3.17 Energy dissipation in ENF tests .................................................................... 65 Figure 3.18 Dimensions of the undulation region (Dimensions in mm) .......................... 66 Figure 3.19 Models with through cohesive layers ............................................................ 67 Figure 3.20 Load-displacement curves of Model D, E and F ........................................... 67 Figure 3.21 Load-displacement curves of Model A (one undulation) and E (two undulations) ............................................................................................................... 68 Figure 4.1 One-dimensional continuum damage mechanics model ................................. 75 Figure 4.2 Three-dimensional progressive damage modeling methodology .................... 79 Figure 4.3 Continuum damage model ............................................................................... 80 Figure 4.4 Elastic softening model with stiffness degradation ......................................... 81 Figure 4.5 Stress-strain curve of the crank band model .................................................... 82 Figure 4.6 Damage evolution in the domain of elastic domain threshold ........................ 84 Figure 4.7 The three damage modes in fiber-reinforced composites ................................ 86 Figure 4.8 UMAT flowchart ............................................................................................. 98 Figure 4.9 Modeling of the quasi-static indentation test in [23] ..................................... 100 Figure 4.10 Load-displacement curve comparison ......................................................... 101 Figure 4.11 Delamination comparison between numerical and experimental results in [23] .......................................................................................................................... 102 Figure 4.12 Delamination with intralaminar damage deactivated .................................. 102 Figure 4.13 Modeling of [Old] laminate .......................................................................... 103 xii Figure 4.14 Model I with the same mesh through thickness ........................................... 104 Figure 4.15 Delamination pattern comparison between experimental (center) from [2] and numerical (left and right) using Model I .................................................................. 106 Figure 4.16 Damage modes in large-angle [0/ a] laminates (u>=30) ............................. 106 Figure 4.17 Model 11 with mesh aligned with ply fiber orientation ................................ 108 Figure 4.18 Result comparison of [0/15] laminate using Model I (left) and Model II....109 Figure 4.19 Damage evolution in [0/15] (Top view) ....................................................... 109 Figure 4.20 Fiber bridging in a [0/10] laminate .............................................................. 1 10 Figure 4.21 The comparison of predicted normalized delamination areas between numerical simulation and mismatch theory [2] ....................................................... 111 Figure 5.1 Organization of dynamic progressive damage analysis ................................. 119 Figure 5.2 VUMAT flowchart ......................................................................................... 122 Figure 5.3 Static tension test on a single element ............................................................ 123 Figure 5.4 Force-displacement comparison ..................................................................... 123 Figure 5.5 Internal energy comparison ............................................................................ 124 Figure 5.6 Transverse impact of alumintun plate model ................................................. 125 Figure 5.7 Comparison of force history with experiment [10] ........................................ 126 Figure 5.8 Comparison of force-displacement curve with experiment [10] .................... 126 Figure 5.9 Comparison of displacement contour with experiment [10] .......................... 127 Figure 5.10 Shell—solid coupling model of [0/90]s plate for static analysis .................... 128 Figure 5.11 Force-displacement curve comparison ........................................................ 130 Figure 5.12 Comparison of damage ................................................................................ 131 Figure 5.13 Damage evolution (only a quarter of specimen shown) .............................. 132 Figure 5.14 Stresses at top 0/90 interface ........................................................................ 133 xiii Figure 5.15 Stresses at bottom 90/0 interface .................................................................. 134 Figure 5.16 Damage sequence in Q3D samples from experiment [13] .......................... 135 Figure 5.17 The generic shell-solid model used for L[O/ 90]6, 2D[O/ 9O]6 , and Q3D[O/ 90]6 specimens ........................................................................................... 135 Figure 5.18 The local solid models ........................................................................... 138 Figure 5.19 The load-deflection curve from experiment [13] ......................................... 143 Figure 5.20 Load-deflection curves from computation ................................................... 144 Figure 5.21 Damage comparison ..................................................................................... 145 Figure 5.22 A quarter of [0/90] laminated plate (a) with damage specified to occur only to the nine elements (b) at the center ........................................................................... 148 Figure 5.23 Post-impact damage of the nine elements (a) using Eqn (5.2) (b) using Eqn (5.3) (0 denotes no damage; damage increases with the value increases) ............... 149 xiv Chapter 1 Introduction 1.1 Background Because of their high stiffness-to-weight ratio, fiber-reinforced polymer-matrix composites have been studied extensively and used successfully in aircraft, ships, high performance automobiles, civil infrastructures, and armor protection [1]. Laminated structures made of unidirectional lamina as shown in Figure 1.1(a) are among the most common forms of fiber-reinforced polymer-matrix composites. Fiber orientation in each lamina and stacking sequence of laminated structures can be chosen to achieve desired strength and stiffness for specific applications. In addition to unidirectional laminates, two-dimensional weaves as illustrated in Figure 1.1(b) are also often laminated together to produce composites with balanced in-plane material properties. However, laminated composites are known for their poor inter-laminar shear strength and poor through-the- thickness properties due to the weak matrix between the plies. This leads to large impact damage, low compression-after-impact strength, and low fracture toughness. To improve the inter—laminar shear resistance, through-the-thickness fibrous reinforcements by means of stitching [2,3], z-pinning [4] and z-binder yarns [5] are introduced into laminated composites. A problem of these reinforcements in the through-the-thickness direction is that they can degrade the in-plane properties [5, 6, 7, 8]. To achieve an optimum balance between in-plane properties and delamination resisance, quasi-three-dimensional (Q3D) woven designs (Figure 1.1(c)) have been developed. Similar to the two-dimensional woven designs, the Q3D woven designs consist of three fundamental constituents: fill and warp tows, and matrix. However, the tows in the Q3D woven designs create interlocking through-the-thickness layer by layer, which eliminates the weak resin-rich interfaces between layers in laminated composites without introducing additional reinforcement in the thickness direction as in the 3D composites. Therefore, it is expected that the Q3D woven composites offer a balanced solution between traditional two- dimensional laminated and 3D fiber-reinforced polymer-matrix composites in terms of in-plane stiffness and delamination resistance. White - Fill yarn Black - Warp yarn Gray - Matrix W ite - 0° lamina Larlninatfd rexglon e— lnterface , e—lnterface - o - r Black 90 lamina Undulation region (a) Laminated (b) Two-dimensional woven White - Fill yarn Black - Warp yarn Gray - Matrix Laminated region 1 f ,I 3‘ \‘ Interfaces between layers eliminated by interlocking I‘m--- . ’t 7‘ s, Undulation region (b) Q3D woven Figure 1.1 Fiber-reinforced polymer-matrix composites studied in this work 1.2 Statement of problem The Q3D woven composites are developed to further improve the delamination resistance of two-dimensional woven composites by eliminating the resin-rich interfaces using through-the-thickness interlocking in the undulation regions. Being different from laminates, both two-dimensional and Q3D woven composites have undulation regions. Furthermore, Q3D composites have slightly larger undulation angles than their two- dimensional counterparts. The objective of this work is to differentiate the Q3D woven designs from the traditional laminated and two-dimensional woven counterparts using numerical models. Given the design differences among Q3D, two-dimensional, and laminated composites, the following questions rise naturally: How is the in-plane stiffness of the Q3D woven designs compared with that of the laminated and two-dimensional woven designs? What are the effects on the delamination propagation introduced by the undulations in woven designs? How does damage initiate and propagate in Q3D composites? Do the Q3D woven designs provide more damage resistance than laminated and two- dimensional woven designs when they are subjected to transverse impact loadings? This work is aimed to answer these questions using finite element methods. 1.3 Research approach and organization To answer the first question above, an average field homogenization method based on finite element analysis of representative unit cells with periodic boundary conditions is presented and implemented in ABAQUS/Standard [9] to predict the effective stiffness of laminated and woven composite designs. The model is validated by experiments. It can then be used to handle complex geometries and perform homogenization of general anisotropic materials. This part of study will be summarized in Chapter 2. To investigate the effects of undulation in woven composites on delamination propagation introduced by the undulations in woven designs, the cohesive zone model (CZM) with plain strain assumptions in ABAQUS/Standard is employed to simulate double cantilever beam (DCB) and end notch flexure (ENF) tests on plain weave models and compare with their corresponding laminated counterparts. Chapter 3 presents the results of the studies. In order to study the damage mechanisms of composites such as the interaction between inter-laminar and intra-laminar damage modes, a three-dimensional Continuum Damage Mechanics (CDM) model for the prediction of the initiation and propagation of intra-laminar damage mechanisms is developed and implemented in ABAQUS/Standard using a user-written material subroutine (UMAT). This three-dimensional CDM model is combined with the Cohesive Zone Model available in ABAQUS/Standard to simulate three-dimensional progressive damage analysis of laminated composite plates subjected to transverse static loadings. Details of the study are given in Chapter 4. The developed three-dimensional Continuum Damage Mechanics (CDM) model is then implemented using a user-written material subroutine (VUMAT), and combined with the Cohesive Zone Model available in ABAQUS/Explicit [10], to compare the dynamic performance of the Q3D woven designs with traditional laminated and two- dimensional woven designs subjected to transverse low-velocity impact. To address the modeling challenges in the progressive damage analysis of woven composites, the shell- solid coupling technique available in ABAQUS is used. This technique offers the numerical accuracy of a full three—dimensional solution and the computational efficiency of shell finite element model. Chapter 5 addresses this part of the study. Chapter 6 summarizes the conclusions from this dissertation research and suggests future studies. REFERENCES . P. Mouritz, M. K. Bannister, P. J. Falzon, and K. H. Leong, Review of applications for advanced three-dimensional fiber textile composites, Composites Part A: Applied Science and Manufacturing, Volume 30, Issue 12, December 1999, Pages 1445-1461 . A. P. Mouritz, K. H. Leong, and I. Herszberg, A review of the effect of stitching on the in-plane mechanical properties of fiber-reinforced polymer composites, Composites Part A 28A (1997) 979-99 . L. C. Dickinson, G. L. Farley, and M. K. Hinders, Prediction of Effective Three- Dimensional Elastic Constants of Translaminar Reinforced Composites, Journal of Composite Materials, 1999; 33; 1002 . A.P. Mouritz, Review of z-pinned composite laminates, Composites Part A: Applied Science and Manufacturing, Volume 38, Issue 12, December 2007, Pages 2383-2397 . P. J. Callus, A. P. Mouritz, M. K. Bannister, and K. H. Leong, Tensile properties and failure mechanisms of 3D woven GRP composites, Composites Part A: Applied Science and Manufacturing, Volume 30, Issue 11, November 1999, Pages 1277-1287 . Y. Yan and YD. Ding Comparison of static properties of 2D and 3D woven fabric reinforced composites. J Compos Technol Res 1999;21:33—6. . L. C. Dickinson, G. L. Farley, and M. K. Hinders, Translaminar reinforced composites: A review, Journal of Composites Technology & Research, JCTRER, Vol. 21, No. l, 1999,pp.3-15 . L. C. Dickinson, G. L. Farley, and M. K. Hinders, Prediction of Effective Three- Dimensional Elastic Constants of Translaminar Reinforced Composites, Journal of Composite Materials, 1999; 33; 1002 . ABAQUS/Standard User’s Manual, online version 6.8 (2008), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. 10. ABAQUS/Explicit User’s Manual, online version 6.8 (2008), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. Chapter 2 Characterization of quasi-three-dimensional woven fabric material 2.1 Abstract This chapter discusses issues related to the characterization of elastic behavior of quasi-three-dimensional (Q3D) woven composite designs. The average field homogenization method based on finite element analysis of representative unit cells with periodic displacement boundary conditions is used at the meso-scale to compute the effective stiffness of Q3D woven composite designs and to compare with the conventional laminated and woven counterparts. Particular attention is given to the effect of the geometric parameters of the unit cell on the in-plane stiffness of the Q3D woven designs. Physical insight is presented for designing Q3D woven composites with high in- plane stiffness. 2.2 Introduction To save computational modeling and analysis efforts of composites with complex geometric microstructures, accurate computation of the effective material properties of composite materials or the so-called homogenization of macroscopic properties is of special interest to researchers in composite material design and characterization. This chapter investigates the elastic behavior of the innovative quasi-three-dimensional (Q3D) woven composite designs and compare their in-plane stiffness with the conventional laminated and woven counterparts. The chapter is organized as follows: First, a literature review of the homogenization methods is presented. The average field homogenization method based on finite element analysis of representative unit cells with periodic boundary conditions is then presented. After this, the implementation scheme in ABAQUS is given. The model is then validated using experimental results from a literature and previous graduate’s work. Parametric studies of Q3D woven composite designs and comparison with conventional laminated and woven counterparts are presented next. Finally, concluding remarks that summarize the results are given in the last section. 2.3 Homogenization literature Review Numerous literature reviews e. g. [1-6] are available on the topic of homogenization of composites. All of the models exploit the periodicity of woven fabric composites to isolate a representative unit cell (RUC). Then the effective material properties are calculated by analyzing the RUC using certain assumptions about the geometry, field distribution, and boundary conditions of the RUC. Earlier attempts at predicting the effective material properties of composites approached this problem from the bounding micromechanical models [7-11]. Depending on the assumptions used c. g. constant state of stress or strain within the composites, some of these models are more accurate than others in giving the lower or upper bounds of the material properties. Though the bounding properties obtained using these models provide a quick estimation of the range of the effective material properties of the composites, they can not be used for accurate FEA of composite structures. Later analytical methods [12-23], e.g. mosaic model, crimp model [15] or sub-cell models [17], utilized classical or modified laminate plate analytical continuum theory methods and simplified geometric modeling of the RUC to compute the effective material properties. Though simple formulas obtained in these methods are good for parametric studies, these analytical methods are either inherently one- or two- dimensional or are only accurate for composites with simple microstructures because of the assumptions used. By dividing the RUC into subcells and calculating the effective material properties for each subcell using analytical methods discussed above, semi-analytical methods [24- 31] eliminate the need to discretize the unit cell in thickness direction; and thus, they give balanced accuracy and efficiency compared with analytical method and full 3D FEA. However, these methods do not provide insight into the non-uniform stress distribution in the constituent materials that characterizes woven composites. Asymptotic Expansion Homogenization [32-38], also frequently called the mathematical homogenization (MH) approach, utilizes a multi-scale perturbation method to establish mathematical relations between the microfields and the macrofields. It is an elegant and systematic approach to examine the micro-mechanics and compute the effective material properties of complex microstructures. Due to a limitation of the homogenization theory, which assumes the unit cell should be periodically represented in any specific area, the through-thickness material properties prediction may be not accurate for thin composites. To overcome this limitation, a modified or enhanced asymptotic homogenization method (EHM)[37,38] is proposed. The average field method [39-45] is based on the fact that the effective mechanical properties measured in experiments are relations between the volume average of the strain and stress of microscopically heterogeneous samples. The average field method [43-45] gives the same effective material properties as the Asymptotic Expansion Homogenization when the representative unit cell is subjected to a series of test strains or test stresses under periodic boundary conditions. The homogeneous displacement and traction boundary conditions provide upper and lower bounds for all possible effective moduli [45] and the periodic boundary conditions provide the effective moduli between these bounds. In [5], the apparent engineering moduli computed using a representative volume element subjected to homogeneous boundary conditions are shown to converge asymptotically to effective moduli of a periodic composite from below and above with increasing number of uniformly—spaced inclusions in the representative volume element. The average field homogenization method based on finite element analysis of representative unit cells with periodic boundary conditions is chosen in this study because of its accuracy. Furthermore, in this study the material properties of the tow and matrix constituents are assumed known and the homogenization is performed at the meso-scale. 2.4 The average field homogenization method 2.4.1 The repetitive unit cell Woven composites can be viewed as periodic structures constructed by the translations of their representative units cells in the three-dimensional domain. For a periodic structure, a representative unit cell can be identified to include all characteristics of the composite such as periodic geometry, material properties, and the periodicity of the field variables; and the periodic structure can be formed by a collection of the translations of the representative unit cell in the directions along which periodicity exists within the solid structure. Herein, the representative unit cell is defined as the “irreducible” unit cell to mean the smallest region capable of building the whole periodic structure only by repeating itself through translations without applying any symmetry and rotation transformations. For the study’s purpose, we choose one such representative unit cell as the reference unit cell, and all other representative unit cells can be obtained by the translation of the reference unit cell. 10 Figure 2.1 illustrates the periodic structures and unit cells of a 2D plain weave (2D) and a 2D 5 harness satin weave (2DH5). The overall response of these periodic composite structures can be characterized in terms of the effective engineering properties, or the so-called homogenized engineering properties, which can be obtained by analyzing their unit cells. (a) Five harness satin weave (2DH5) structure (b) Unit cell of 2DH5 Figure 2.1 Composite structures and their unit cells 2.4.2 The average strain theorem for periodic boundary conditions After the representative unit cells of the composite structures are selected, the appropriate boundary conditions must be applied to the unit cells in order to produce accurate effective engineering properties of the composite structures. The most often used boundary conditions in literature include the homogeneous displacement boundary 11 conditions, homogeneous traction boundary conditions, and periodic boundary conditions. Discussions about the selection of appropriate boundary conditions can be found in [5, 39, 44]. Periodic boundary conditions are chosen for this study because of their accuracy. It has been shown [11, 44] that homogeneous boundary conditions applied on the surface of a homogeneous body will produce a homogeneous field. Specifically, the homogeneous displacement boundary conditions associated with constant test strains applied on the surface of a representative unit cell will produce average strains in the representative unit cell identical to the applied constant test strains. This is known as the average strain theorem. Similarly, the average stress theorem states that the homogeneous traction boundary conditions associated with constant test stresses applied on the surface of a representative unit cell will produce average stress fields in the representative unit cell identical to the applied constant test stresses. Here, the average strain theorem is shown to also hold for periodic displacement boundary conditions. Using this theorem, efforts can be saved to calculate the average strains by avoiding volume integration over the whole unit cell. In the following, the constant test strains 83- will be shown to be identical to the average strain EU for periodic displacement boundary conditions. For a two-phased composite with tow material denoted by “1” with volume V] and bounding surface S1, and matrix material denoted by “2” with volume V2 bounding surface 52. SI and S2 consist of the interface S12 between the tow phases and the external surface S. Substituting 1 8ij =§(ui,j+uj,i) (2.1) 12 into - 1 817 = — [5‘ng (2.2) V v yields 2vK = [V1 (uf‘). + u“). )dV + [V2 (1.4%) + 1.33,th (2.3) Then applying Gauss theorem IV uLJ-dV =15 uinde (2.4) it becomes that 2vE,-,- = [Sr (1.9),: + u(l)n,-)dS + 152(ujzln + u(2)n )dS (2.5) Assuming perfect contact between the phases with interfaces S12 , i.e. up) =ufz) on 512 (2.6) it leads to the cancellation of the contributions from 512. Therefore, we have - 1 8;} =st(uinj +ujni)dS (2.7) The external surface S can be divided into boundary pairs with coordinates X0 and X0 +l , where ldenotes the translation/periodicity vector of the unit cell. Applying u (X0 +l) —u (X0): 8,)- lj (2.8) to the opposite boundary surface pairs of the unit cell gives 2.,- = 83 (2.9) 13 This shows that the average strain theorem holds for periodic displacement boundary conditions. Using the displacement and material periodicity, strain-displacement relationships, and material constitutive laws, the periodic displacement boundary conditions can be shown to be sufficient to guarantee periodic strain and stress distribution for periodic materials in displacement-based finite element formulations. In a periodic model, it is assumed that periodicity extends to infinity in all the directions along which periodicity exists within the solid structure. The correct boundary conditions on a unit cell should lead to a periodic distribution of the field quantities, namely, displacement, strain, and stress across the whole periodic structure. The periodic conditions require that the displacements in the various unit cells differ only by a constant, and the strains and stresses are identical in all of the unit cells. The relationship of displacement, strain and stress components between a master node and its image nodes in a unit cell can be expressed as follows [39]: EU (X0 +1) =8U(X0) (2.10) al-J-(Xo +1) =a,-j(xo) (2.11) u,(X0 +1) =u,(xo) +2.71]. (2.12) where lj is the periodicity vector component linking a master node to its image nodes in a unit cell. Using periodicity and constitutive laws, it is shown below that periodic strain fields (Eq. (2.10)), and stress fields (Eq. (2.11)) are guaranteed by specifying the periodic displacement fields (Eq. (2.12)) for periodic structures. Assuming a prescribed array of 14 constant test strains 83 applied to the corresponding nodes on the opposite boundaries of the unit cell, which are identical to the average strain E17 by the average strain theorem (Eq. (2.9)), then Eq. (2.11) gives 91 ul-(X0+l) =ui(X0)+81jj (2.13) . 1 . . . Since 3,-1- = E(u i, j+u J'J ) , and £8 and lj are constant for a given periodlc structure and displacements are continuous, we have 80- (X0 + l) =£ij(X0) (2.14) Because of the periodicity of material property, the following holds: Cijk/(XO +1) = Cijkl (X0) (2.15) Applying the stress-strain relationship: 0,-1- = C ijkl 8k, yields 0ij(X0 + l) =0'iJ(X0) (2.16) Therefore, the periodic displacement fields (Eq. (2.12)) for periodic structures yields periodic strain (Eq. (2.10)) and stress fields (Eq. (2.11)). For a general anisotropic material, the 6x6 stiffness matrix and the corresponding compliance matrix have 36 entries. To find the 36 entries, periodic displacement boundary conditions corresponding to six independent constant test strains are applied to the boundaries of the unit cell. These six test strains will yield six average stresses. Thus the 36 entries can be computed using the constitutive relationship a = Ea (2.17) 15 where E is the average stress associated with the average strain E , and C is the homogenized stiffness matrix. Using the six independent constant test strains given below: 80(1) = {1,o,o,o,o,0}T 30(2) = {0,1,0,0,0,0}T 30(3) = {0,0,1,0,0,0}T 80(4) = {0,0,0,1,0,0}T 30(5) = {0,0,0,0,1,0)T 20(6) = {0,0,0,0,0,1}T (2.18-2.23) and the average strain theorem, the homogenized stiffness matrix can be constructed conveniently as follows: C : [3(1),E(2),E(3),5(4),5(5),Em)] (2.24) where 6(i),i=1,2,...6 are the averaged stresses corresponding to the six independent constant test strains so“) given by Eqs. (2.18-2.23). Explicitly, the stiffness matrix can be expressed as follows: PC” C12 C13 C14 C15 El6q 01d) 01(12) 01(13) 6(4) (5) (6) C21 C22 C23 C24 C25 C26 0'22 322 3 — _ _ — — - -1 -2 -3 -4 -5 C31 C32 C33 C34 C35 C36 03(3) 0353) 03(3) 03(3) 03(3) 333 C41 C42 C43 C44 C45 C46 31%) 53) 51(3) 31(3) 31(3) 51(3) 0| 11 55] C52 C53 C54 C55 C56 31%) 3632) 51(3) 51g) 31(35) 31(3) _C61 562 563 C64 C65 C66] —(1) 3(2) —(3) 3(4) —(5) 3(6) 16 By the inversion of the stiffness matrix yields the compliance matrix, the homogenized compliance matrix is obtained as follows: - ——1 S = C (2.26) Strictly speaking, woven composite materials are not orthotropic because they do not generally possess three mutually perpendicular planes of material symmetry. But they can be approximated as orthotropic and characterized by nine independent effective elastic constants. The nine elastic constants (E11,E22,E33,512,Gl3,523,512,F13,523) for an orthotropic material are obtained by the inversion of the following relationship: F l/Ell —I/-21/-E-22 -F31/E33 0 0 O Jig/E11 1/E22 532/1533 0 0 0 §= —;13/E11 —1/-23/E22 l/E33 0 0 O o o 0 1/012 0 0 0 o 0 0 1/513 0 _ 0 o 0 o 0 1/523 (2.27) 2.5 Implementation in ABAQUS/Standard Using displacement-based finite element formulations in ABAQUS/Standard [46], the periodic displacement boundary conditions (Eq. (2.12)) can be applied to the corresponding nodes at the boundaries of any parallelepiped unit cells using keyword *EQUATION. The flowchart in (Figure 2.2) illustrates the implementation procedure used in this study. In order to implement the homogenization procedure described above, several programs have been developed. The procedure is given as follows: A PYTHON Pre- Processing script can be called using the “abaqus cae script” command to model and mesh the parts and generate ABAQUS input files automatically. Then, a MATLAB Pre— 17 Processing program is called to generate node pairs on the unit cell for applying periodic boundary conditions. After this, the ABAQUS input files are modified manually to apply the periodic boundary conditions. Then a general static analysis is performed using the ABAQUS/Standard solver to generate output database files, which are accessed by calling a PYTHON Post-Processing script to extract the desired filed variables. Finally, the MATLAB Post-Processing program is used to compute the effective stiffness. PYTHON Pre-Processing Script Create model, generate ABAQUS input files I MATLAB Pre-Processing Program Generate node pairs on unit cell boundaries I Modify ABAQUS input files to apply periodic boundary conditions I Run ABAQUS Analysis Generate output databases PYTHON Post-Processing Script Generate result files (i.e. stress, strain) from ABAQUS output databases 7 I MATLAB Post-Processing Program I Perform homogenization to compute ! effective stiffness Figure 2.2 Implementation procedure of the homogenization method 18 2.6 Model calibration using a plain weave problem To calibrate the developed model, the nine effective engineering constants of a plain weave problem computed using this model with periodic displacement boundary conditions are compared with the experimental and analytical results in [30]. (a) Meshed unit cell 2h 2a w (b) Unit cell dimensions for two harness woven Figure 2.3 Unit cell of plain weave The unit cell used is illustrated in Figure 2.3, and its dimensions are summarized in Table 2.1. The mesh seeds used are as follows: 6 elements for length a0 , 4 elements for length a“ , and 2 elements for ply thickness h. The undulation length of the unit cell au is calculated based on the linear geometry and the fiber volume fraction provided. The material properties for tow and matrix are given in Table 2.2 and Table 2.3. As shown in Table 2.4, the results obtained using the method described above correlate well with the experimental results except for the in-plane shear modulus. The deviation in the effective 19 in-plane shear modulus may be explained by the difference of the boundary conditions imposed between the numerical model and experiment. As discussed in [44], restricting deformed unit cells of composites to remain a parallelogram with straight edges for shear loading is an overly restrictive constraint, and it yields larger effective stiffness values. Table 2.1 Plain weave unit cell geometry [22] Half unit cell Tow width Undulation Ply thickness Unit cell length 510 (mm) length h (mm) thickness a (m) au (mm) 2h (mm) 0.87 0.60 0.27 0.05 0.10 Table 2.2 Tow Material Properties [22] Material Err E22 =5 33 012 = G13 G23 V12 = V13 V23 (GPa) (GPa) (GPa) (GPa) E-glass/Vinylester 57.50 18.80 7.44 7.26 0.25 0.29 vf = 0.80 Table 2.3 Matrix Material Properties [22] | Material l E (GPa) | G (GPa) v | IVinylester _| 3.40 | 1.49 0.35 J Table 2.4 Comparison of homogenized material properties E11 =E22 E33 512 513 = 523 F12 F13 = 523 (GPa) (GPa) (GPa) (GPa) Experiment 24.8 11.1 85:26 6.5208 4210.7 1 011001 0281007 [221 Analytical 25.33 13.46 5.19 5.24 0.12 0.29 122] Current 25.13 12.78 4.35 3.88 0.15 0.30 20 2.7 More validation with experiments In this section, the developed model is further validated using cross-ply laminates L[0/90]6, 2D plain weave 2D[0/90]6, two-hamess [0/90]6, and five harness Q3DH5[0/90]6 from the experimental work [47]. 2.7.1 The specimens The unit cells (see Figure 2.4) are modeled with 12 plies through-the-thickness which are the same as in the experimental specimens. The dimensions of the unit cell are given in Table 2.5. Because the unit cells use the same thickness as in the experimental specimens, periodic boundary conditions are applied only in the in-plane directions of the unit cells and homogeneous boundary conditions are applied on the top and bottom surfaces. This type of boundary condition is denoted as IPBC. To study size effect, multiple cell specimens (see Figure 2.5) matching the dimensions of the experimental specimens in three dimensions are analyzed with homogeneous boundary conditions in three directions. This type of boundary condition, (a) L[O/90]6 (b) 2D[0/90]6 (c) Q3D[0/90]6 (d) Q3DH5 [O/9O]6 is denoted as HBC. Figure 2.4 Unit cells analyzed with 12 plies 21 Table 2.5 Unit cell geometry [47] Half unit cell Tow width Undulation length Ply thickness length a0 (mm) a“ (mm) h (m) a (mun) 9.96 9.30 . 0.66 0.17 (a)L[0/90]6 (4 x4 cells) (b) 2D[0/90]6 (4x4 cells) 0 (c) Q3D[O/90]6 (4x4 cells) ((1) Q3DH5[O/9O]6 (2x2 cells) Figure 2.5 Specimens analyzed to match the dimensions in experiments 2.7.2 Material properties The experimental specimens are made of CYCOM 1003/W-490 E-Glass Prepregs with resin weight fraction 36% (equivalent to 45% fiber volume fraction) and matrix from Cytec Industries Inc. To make thin specimens, the woven specimens are cured and processed. Therefore, the fiber volume fraction and material properties are changed. Due to the lack of the material properties of tows and matrix of the specimens in [47], moduli of the CYCOM 1003/W-490 E-Glass Prepregs and matrix from Cytec Industries Inc. (see Table 2.6) are scaled by the ratio (13.80/(38.61+8.27) ) of the cross ply laminate stiffness between the reported result in [47] and the calculated result using the above material properties to obtain the input tow material properties as shown in Table 2.8. 22 Table 2.6 Tow material properties CYCOM 1003/W-490 E-glass Material I 511 J 1522:1333 G12=Gi3 G23 V12=V13 V23 (GPa) (GPa) (GPa) (GPa) IE-lass/EpoxyL38.6l | 8.27 | 4.14 | 3.43 | 0.26 |0.40| Table 2.7 Matrix material properties from Cytec Industries inc. Material | E (GPa) | v ] Epoxy | 3.45 | 0.35 | Table 2.8 Tow material properties used in the analysis (Interpolated) Material Err 522 =1533 G12 = 013 G23 V12 = V13 V23 (GPa) (GPa) (GPa) (GPa) | E—lass/Epoxy | 22.57 L 4.84 | 2.42 J 2.01 | 0.26 ] 0.40 | 2.7.3 Results and discussions The computed effective material properties are compared with the experimental results from [47] as given in Table 2.9. The experimental results are from compression test and only in—plane stiffness is reported. The homogenized material properties show the same trend as the experimental results, i.e. the stiffness of the composites are in the following order: L>2D>Q3DH5>Q3D Apparently, the laminated composite has the largest in-plane stiffness because it has neither matrix pockets nor undulation in it. The in-plane stiffness of the Q3DH5 weave is always greater than its Q3D counterpart because the Q3DH5 weave has less undulation than its Q3D counterpart. As it is shown later in the parametric studies, Q3DH5 weaves can have larger in-plane stiffness than their 2D plain weave counterparts because stiffness changes with the change of the undulation angle and fiber volume fraction. 23 The computed effective material properties for the woven composites are consistently about 5% lower than the experimental results. This might be due to the interpolated tow material properties used. In this case, the results show no improvement in the homogenized results using the actual size of the experimental specimen than using only one unit cell. Table 2.9 Comparison of results with experiments E11 =E 22 (GPa) Tested specimen Experiment Numerical Error Numerical [47] IPBC (unit cell) (%) HBC(Actual size) L [0/90]6 13.80 13.79 -0.07 13.79 (4x4 cells) 2D[0/90]6 13.63 11% 12.98 -4.80 12.97(4x4 cells) Q3D[0/90]6 13.31 :t6% 12.69 -4.68 12.68(4x4 cells) Q3DH5[0/90]6 13.62 12.95 -4.93 12.94(2x2 cells) 2.8 Parametric study This section investigates the effects of the unit cell dimensions on the effective in- plane stiffness of the 2D and Q3D woven designs by performing parametric studies using the developed method. 2.8.1 The unit cells and material properties used For the parametric studies, the unit cells analyzed for different woven designs are illustrated in Figure 2.6. The dimension diagram for two harness woven designs is the same as illustrated in Figure 2.3 (b), and the dimension diagram for five harness woven designs is given in Figure 2.7. Unit cells with two plies in thickness direction are used for the 2D weaves, while unit cells with six plies ( the minimum number of plies in the Q3D weaves) are used for Q3D and Q3DH5. Periodic boundary conditions are imposed in three directions of the unit cells. The mesh seeds used are the same as before: 6 elements for length a0, 4 elements for length au, and 2 elements for ply thickness h. 24 The input material properties used are given in Table 2.10 and 2.11. The effective moduli are normalized with respect to the effective moduli of [0/90]s laminate given in 9 (a) 2D [0/90] (b) 2DH5 [0/90] Table 2.12 . (c) Q3D[0/90]3 (d) Q3DH5[0/90]3 Figure 2.6 Unit cells used for parametric study -- "1" ------ 6h 1 1 1r . ....... a . .2. 4.13.. L 5a Figure 2.7 Unit cell for five harness woven 7 Table 2.10 Tow 7Material Properties [43] 7 - 1 L Mammal I E11 1 E22 =E33 I G12 =G13 G23 v12 =V13 V23 (GPa) (GPa) (GPa) (GPa) firaghite/Epoxyl 134.00] 10.20 | 5.52 [3.43] 0.30 |0.49| 7 Table 2.11 Matrix Material Properties [43] 7 [ Material L E (GPa) | G (GPa) | v | | Epoxy | 3.45 | 1.28 | 0.35 | 7 Table 2.12 [0/90]s laminate effective moduli [43] 7 LMaterial I 1511:?22 (GPa) I E33 (GPa) I G12=G13(GPa) I "623 (GPa) I [0/90]s 72.50 10.20 5.52 4.26 laminate 7 25 2.8.2 Parametric study by varying a/h For the parametric studies by varying the ratio of the half cell length over ply thickness “a/ h ” (see Figure 2.7), the dimensions of the unit cell are summarized in Table 2.13. h is scaled according to the ratios a/h =10, 25, 50 and 100, while a is kept constant. Table 2.13 Unit cell dimensions for varying a/ h Half unit cell Tow width Undulation length Ply thickness length a0 (mm) au (mm) 1'1 (mm) a (m) 10 8.33 1.67 0.01~1.00 The normalized in-plane effective moduli denoted as E1 1 = E22 and G 12 by varying a/ h are plotted in Figure 2.8. Based on the results, the following conclusions can be made: 0 Of all the designs, in-plane stiffness increases with the increase of a/ h because the undulation angle decreases with the increase of a/ h . 0 Higher harness woven designs have higher in-plane Young’s Moduli due to the lower number of undulations in higher harness woven designs. 0 2D weaves have higher in-plane stiffness than their corresponding Q3D counterparts due to smaller undulation angles in 2D weaves. The differences become smaller with the increase of “ a/h 0 The engineering constants of all woven designs do not approach those of laminated counterparts because of the existence of the pure matrix zones and the undulation regions in the woven designs 26 _ ,- - -o- - 20 0' .’ —o—201-15 ' - -A- - 030 .. “Aswww .__._..__3-L.__._,,_.__.._ _ L_ _ 0 40 +Q3DH5 0.20 ———-e *4 4 ~ gm, ,3 we “—- m~ 0.00 r r . r . 0 20 4o 60 80 100 - «O- - ZD —O——ZDH5 - -A- - Q3D +Q3DH5 (b) Shear Modulus Figure 2.8 Normalized in-plane effective Moduli by varying a/h 27 REFERENCES . P.W. Chung and K.K. Tamma, Woven fabric composites — developments in engineering bounds, homogenization and applications. Int J Numer Meth Engng 45 (1999), pp. 1757—1790. . J. Bystrom, N. Jekabsons, and J. Varna, “An evaluation of different models for prediction of elastic properties of woven composites”, Composites: Part B 31 (2000) 7—20 . A. E. Bogdanovich, “Multi-scale modeling, stress and failure analyses of 3-D woven composites” J Mater Sci (2006) 41:6547—6590 L. Leon, Mishnaevsky Jr and Siegfried Schmauder “Continuum mesomechanical finite element modeling in materials development: A state-of-the-art review” Appl Mech Rev vol. 54, no 1, January 2001. A. Drago and M.J. Pindera,"Micro-macromechanical analysis of heterogeneousmaterials: Macroscopically homogeneous vs periodic microstructures”, Composites Science and Technology 67 (2007) 1243 - 1263 . L. Onal and S. Adanur, Modeling of Elastic, Thermal and Strength/Failure Analysis of Two-Dimensional Woven Composites—A Review, Applied Mechanics Reviews, January 2007, Vol. 60/37 . Z. Hashin, The elastic moduli of heterogeneous materials. Journal of Applied Mechanics Transactions of the ASME 1962; 29E: 143-150. . Z. Hashin and S. Shtrikman A variational approach to the theory of the elastic behaviour of multiphase materials. Journal of the Mechanics and Physics of Solids 1963; 11:127-140. . Z. Hashin and BW. Rosen, The elastic moduli of fiber reinforced materials. Journal of Applied Mechanics 1964; 31:223 -232. 10. R. Hill, Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids 1963; 11:357-372. 11. J. Aboudi, Mechanics of composite materials: a unified micromechanical approach, Amsterdam, New York : Elsevier, 1991 12. EN. Cox and G. Flanagan, Handbook of analytical methods for textile composites, NASA CR 4570, NASA Langley Research Center, Virginia, USA, 1997 13.8. Z. Sheng and S. V. Hoa, Three Dimensional Micro-Mechanical Modeling of Woven Fabric Composites, Journal of Composite Materials, Vol. 35, No. 19/2001 28 14. J. Whitcomb and X. Tang, “Effective Moduli of Woven Composites, Journal of Composite Materials”, Vol. 35, No. 23, 2127-2144 (2001) 15. T. Ishikawa and T.W. Chou. Thermoelastic analysis of hybrid fabric composites. Journal of Materials Science 1983; 18: 2260 -2268. 16. T. Ishikawa, M Matsushima, Y Hayashi, and T.W. Chou, Experimental confirmation of the theory of elastic moduli of fabric composites. Journal of Composite Materials 1985; 19:443-458. 17. T.W. Chou and T. Ishikawa, Analysis and modeling of Two-Dimensional Fabric Composites, Composite Material Series 3, Textile Structure Composites, T.W. Chou and F.K.Ko, Eds, Elsevier Science Publishers B.V., 1989, pp. 210-264 18. 1.5. Raju and J.T. Wang, Classical laminate theory models for woven fabric composites. Journal of Composites Technology and Research 1994; 16(4):289-303. 19. N.K. Naik and RS. Shembekar Elastic behavior of woven fabric composites: I | Lamina Analysis. Journal of Composite Materials 1992; 26(15):2196 -2225. 20. N.K. Naik and V.K. Ganesh An analytical method for plain weave fabric composites. Composites 1995; 26:281-89. 21. B.V. Sankar and RV. Marrey Analytical method for micromechanics of textile composites. Composite Science and Technology 1997; 57:703 -713. 22. D. Scida, Z. Aboura, M.L. Benzeggagha, and E. Bocherens, A micromechanics model for 3D elasticity and failure of woven-fibre composite materials, Composites Science and Technology 59 (1999) 505—517. 23. S. C. Klann, Unit cell analysis of composite microstructures, M.S. thesis 2008, Michigan State University. 24.13. A. Bednarcyk, “Discussion of “Woven fabric composite material model with material non-linearity for non-linear finite element simulation” by Tabiei and J iang [Int. J. Solids Struct. 36 (1999) 2757—2771]”, International Journal of Solids and Structures ,Volume 38, Issues 46-47, November 2001, Pages 8585-8588 25. A. Tabiei and Y. Jiang, Woven fabric composite material model with material nonlinearity for nonlinear finite element simulation, International Journal of Solids and Structures 36 (1999) 2757-2771 26. B. A. Bednarcyk and M-J. Pindera Inelastic Response of a Woven Carbon/Copper Composite— Part II: Micromechanics Model 27. A. Gonzalez, E. Graciani, and Federico Paris, Prediction of in-plane stiffness properties of non-crimp fabric laminates by means of 3D finite element analysis, Composites Science and Technology 68 (2008) 121 - 131 29 28. V.R. Aitharaju and RC. Averill, Three-dimensional properties of woven-fabric composites, Composites Science and Technology 59 (1999) 1901-1911 29. A. Tabiei and W. Yi, Comparative study of predictive methods for woven fabric composite elastic properties, Composite Structures 58 (2002) 149 - 164 30. R Tanov and A. Tabiei, Computationally efficient micro-mechanical woven fabric constitutive, J Models Appl Mech 2001, 68(4):553 - 60. 31. N. Abolfathi, A. Naik, G. Karami, and C. Ulven, A micromechanical characterization of angular bidirectional fibrous composites Computational Materials Science, 2008 43 (4): 1193-1206 32. P.W. Chung, K.K Tamma, and RR Namburu, Asymptotic Expansion Homogenization for Heterogeneous Media: Computational Issues and Applications, Journal of Composites: Part A — Applications, 2001, 32 (9): 1291-1301. 33. J .M. Guedes, and N. Kikuchi, Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Methods Appl Mech Eng 1990, 83:143-98. 34. N. Takano, Y. Uetsuji, Y. Kashiwagi, and M. Zako, “Hierarchical modelling of textile composite materials and structures by the homogenization method”. Modelling Simul. Mater. Sci. Eng. 7 (1999) 207-231. 35. Y. Ohnishi, M. Zako, and K. Nishiyabu, “Microstructure-based deep-drawing simulation of knitted fabric reinforced thermoplastics by homogenization theory”. International Journal of Solids and Structures,Volume 38, Issues 36-37, 2001, Pages 6333-6356 36. N. Takano, M. Zako, and Y. Ohnishi, Macro-micro uncoupled analysis procedure for nonlinear structural mechanics of composites - RVE approach to microstructural geometrically nonlinear problem, Japan Society of Materials Science, Journal. Vol. 45, no. 2, pp. 163-168. Feb. 1996 37. F. Rostam-Abadia, C-M Chen, and N. Kikuchi, “ Design analysis of composite laminate structures for light-weight armored vehicle by homogenization method”, Computers and Structures 76 (2000) 319-335 38. C-M Chen, N. Kikuchi, and F. Rostam-Abadi, “An enhanced asymptotic homogenization method of the static and dynamics of elastic composite laminates” Computers and Structures 82 (2004) 373-3 82 39. J. D. Whitcomb, C. D. Chapman, and X. Tang, Derivation of Boundary Conditions for Micromechanics Analyses of Plain and Satin Weave Composites, Journal of Composite Materials, 2000; 34; 724 30 40. X. Tang and J. D. Whitcomb, General Techniques for Exploiting Periodicity and 41. 42. 43. 45. 46. 47. Symmetries in Micromechanics Analysis of Textile Composites, Journal of Composite Materials, Vol. 37, No. 13/2003 Z. Xia, Y. Zhang, and F. Ellyin, A unified periodical boundary conditions for representative volume elements of composites and applications International Journal of Solids and Structures, Volume 40, Issue 8, April 2003, Pages 1907-1921 Z. Xia, C. Zhou, Q. Yong, and X. Wang, On selection of repeated unit cell model and application of unified periodic boundary conditions in micro-mechanical analysis of composites, International Journal of Solids and Structures, Volume 43, Issue 2, January 2006, Pages 266-278 Whitcomb JD. Three dimensional stress analysis of plain weave composites. Paper presented at the 3rd Symposium on CompositeMaterials: Fatigue and Fracture, Orlando (FL), 1989. p. 417—38. . C. T. Sun and R. S. Vaidya, Prediction of composite properties from a representative volume element, Composites Science and Technology 56 (1996) 171-179 M. Hori and S. Nemat-Nasser, On two micromechanics theories for determining micro-macro relations in heterogeneous solids, Mechanics of Materials 31 (1999.), 667—682. ABAQUS/Standard User’s Manual, online version 6.8 (2008), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. K. Rosario, Impact Properties of Quasi-3D Composites, M.S. thesis 2008, Michigan State University. 31 Chapter 3 Two-dimensional progressive delamination analysis of plain woven composites 3.1 Abstract This chapter investigates the undulation effects on the interlaminar damage (i.e. delamination) resistance in woven composites using two-dimensional Cohesive Zone Model (CZM) with plain strain assumptions. Numerical simulation of DCB tests on both symmetric and unsymmetric plain weave models and their corresponding laminated counterparts are performed, the test results are compared, and the undulation effect in plain weave models are discussed. To study the effects of undulation direction and number on the interlaminar damage in woven composites, ENF tests are performed on symmetric plain weave models with cohesive layers inserted at different locations; and the test results are compared and discussed. 3.2 Introduction This chapter is aimed to study the undulation effects on the interlaminar damage (i.e., delamination) resistance in woven composites using two-dimensional plain weave models. Progressive delamination will be investigated by simulating DCB and ENF tests using the two-dimensional Cohesive Zone Model (CZM) in ABAQUS/Standard [l] with mixed-mode criteria for delamination initiation and propagation. Special attention is given to the effects on the delamination propagation of undulation in woven composites. In this study, interlaminar delamination is considered. Each ply is modeled as elastic; and thus, intra-laminar damage modes such as matrix cracking are not considered. Furthermore, two-dimensional plain strain is assumed. Thus, the delamination involves at most Mode I and Mode II. 32 3.3 Literature review 3.3.1 Literature review on delamination Fiber-reinforced polymer composites (FRPCs) can accumulate damage before structural collapse [2]. Dominant damage modes in FRPCs include intralaminar damage, e.g., fiber breakage, matrix cracking, matrix crushing, and interlaminar damage, i.e., delamination, which is the separation between adjoining layers. Delamination is one of the most dangerous damage modes in FRPCs, and the prediction of delamination in FRPCs is challenging because of the anisotropy and heterogeneity natures of FRPCs. Reviews of the testing and analysis methodologies of delamination can be found in [3-5]. Because of the material anisotropy and geometric complexity of FRPCs, delamination in FRPCs may involve three modes of fracture simultaneously. The three fracture modes are illustrated in Figure 3.1 (a), (b), and (c) [5]: the interlaminar tension (Mode 1), shearing (Mode II), and tearing (Mode III). Each mode is associated with a fracture toughness value. To characterize the fracture toughness for Mode I delamination, the double cantilever beam (DCB) test as shown in Figure 3.1 (d) is used. The three-point bending end notch flexure (ENF) as shown in Figure 3.1 (e) and end loaded split (ELS) as shown in Figure 3.1 (f) are often used for Mode 11. For two-dimensional studies, the characterization method for Mode III fracture toughness is omitted. In Figure 3.1 (d), (e) and (f), a denotes the beam length, 00 denotes the initial crack length, h designates the height of half a beam, and F denotes the load. 33 m (a) Interlaminar tension (Mode 1) (b) Shearing (Mode 11) (c) Tearing Mode IlI) F do A 2h a F (d) DCB test 610 _I L F a/2 VI A 2h a A I‘ 'l (e) ENF test Rigid clamped end F A a0 L211 A a ‘ (f) ELS test Figure 3.1 The fracture modes and test methods 34 3.3.1.1 Fiber orientation effects on delamination resistance in laminate composites The application scope of the ASTM E399 standard Mode I DCB test method is currently limited to delamination propagation through mid-plane of unidirectional specimens with embedded delaminations [3], and an international consensus of the Mode II ENF test method has not been reached. The validity of the fracture toughness values obtained from these unidirectional specimens with embedded delaminations is debatable because delamination may never propagate between plies of the same orientations due to the lack of stiffness mismatch at the interfaces [6]. Therefore, multi-directionally laminated specimens have been investigated intensively [7-10]. The conclusions of the effects of fiber orientation on delamination from these studies are quite contradictory. Using sub-laminate specimens, Tao and Sun [10] found that the Mode II interlaminar toughness of 0/0 interfaces in AS4/3501-6 composite material decreases as the off-axis angle 0 increases. The test results in [9] showed that the interlaminar fracture energies in unidirectional laminates were significantly lower than those in multidirectional laminates under Mode I, Mode 11, and Mixed-Mode I/II loading conditions. However, Morais et al. [11] observed no significant dependence of the Mode 1 delamination resistance on the fiber orientation. 3.3.1.2 Fiber orientation effects on delamination resistance in woven composites In addition to the study of delamination in laminated FRPCs, experimental investigation of interlaminar delamination in 3D orthogonal interlocked fabric composites [12,13] has also been reported. They found that the Mode I fracture toughness values of composites reinforced in the through-thickness direction by weaving 35 or stitching were higher than traditional composites by factors of nearly 2 and 3 respectively. Chen et al. [15] investigated the delamination in braided T-piece composite specimens using the cohesive zone model in ABAQUS and found that the failure in T- piece composites is dominated by the propagation of the crack rather than its initiation. However, to the authors’ best knowledge numerical studies of the undulation effects on delamination propagation in woven designs have not been reported. The aim of this work is to study the undulation effects on the delamination resistance in woven composites using two-dimensional plain weave models. 3.3.2 Literature review on progressive damage analysis Different approaches have been proposed to model the progressive damage analysis of FRPCs. The multi-level and meso-finite element modeling approaches [15-18] were shown to make good predictions of the damage initiation using strength based failure criteria but make poor predictions of damage propagation. In these models, by assuming periodicity, in which a representative unit cell can be identified to include all characteristics of the composites such as periodic geometry and material properties and uniform damage distribution across whole structures; damages characterized by the degradation of material properties were considered sequentially on representative unit cells from micro scale (fiber/matrix), meso-scale (tow, matrix pocket), and macro-scale (composite structures). However, once material damage is localized into a macrocrack in a narrow band, the assumptions of periodicity and uniform damage distribution are not valid anymore. This leads to the poor predictions of damage propagation in these methods. 36 Fracture mechanics based Virtual Crack Closure Technique (VCCT) [19,20], originally proposed by Rybicki and Kanninen[21], assumes that the energy released during delamination propagation equals the work required to close the crack back to its original position. Based on this assumption, the single-mode components of the energy release rate are computed from the nodal forces and nodal relative displacements. Delamination growth is predicted when the delamination criterion in terms of the critical energy release rates is met. The VCCT approach is computationally effective since the energy release rates can be obtained from only one analysis. However, this approach relies on the assumption of existing initial defects or cracks. Furthermore, it requires complex mesh techniques to advance the crack front when the local energy release rates reach a critical value. The cohesive zone models (CZM) [22-28], originally proposed in [29], are particularly advantageous for delamination modeling of laminated composites where the interfacial strength is relatively low compared with that of the in-ply materials [30,31]. The CZM approach relates traction to separations at an interface where a crack may occur. Damage initiation is related to the interfacial strength, i.e., the maximum traction on the traction-separation curve. When the area under the traction-separation curve is equal to the fracture toughness, the traction is reduced to zero and new crack surfaces are formed [27]. The CZM approach combines the advantages of the strength based failure criteria for the prediction of damage onset and the fracture mechanics based criteria. for crack propagation after damage onset. Thus, it has the capability to predict both onset and propagation of delamination without previous knowledge of the crack location and propagation direction. The CZM has an advantage over stress-based methods because it 37 uses traction-separation laws to overcome spurious mesh dependence. It also has an advantage over the fracture mechanics based approach like VCCT because it does not rely on the assumption of existing initial cracks. Therefore, CZM is widely used to simulate the resin-rich layers connecting adjacent laminae of laminate structures for the progressive delamination analysis. This work investigates the undulation effects on delamination propagation in plain woven designs using the CZM in ABAQUS/Standard. In this study, only delamination is considered. Each ply is modeled as elastic; and thus, intralaminar damage modes are not considered. Furthermore, two-dimensional plain strain is assumed. The out-of-plane shear stress associated with Mode III is assumed zero. Thus, the delamination involves at most Mode I and Mode 11. Progressive delamination will be investigated by simulating DCB and ENF tests using the two-dimensional CZM with mixed-mode criteria for delamination initiation and propagation. Special attention is given to the effects on the delamination propagation of undulation in woven composites. 3.4 Delamination modeling Delamination occurs between adjacent layers of laminated composites. An interface can be inserted between any two adjacent layers where delamination may occur [32]. Using CZM, cohesive elements are inserted between layers of solid elements where delamination is expected to occur. Cohesive elements enable the modeling of fracture initiation and propagation in finite element analyses. However, to obtain a successful finite element simulation using CZM, the cohesive zone contribution to the global compliance should be small enough to avoid the introduction of artificial compliance to the model. Also, the progressive damage analysis using CZM poses numerical difficulties 38 associated with the proper definition of the stiffness of the cohesive layer, the requirement of extremely refined meshes, and the convergence problems caused by material softening. The material parameter determination and numerical issues related to the implementation of CZM are discussed in this section. 3.4.1 Damage evolution law CZM uses damage evolution laws, which are traction-separation curves, to control the crack initiation and growth. Most of the available traction-separation laws assume initially linear elastic traction-separation displacement behavior followed by the initiation and evolution of damage. The most often used traction-separation curves can be found in [27]. The traction-separation relationship for fracture simulation using cohesive finite elements is characterized by fracture toughness, initial stiffness, and strength of the interface. Figure 3.2 illustrates the bilinear damage evolution laws under pure mode conditions [3 1]. In Figure 3.2, kl- and If (i=1, 2, 3) are the interlaminar normal and shear stiffnesses and strengths, di(i=1, 2, 3) are the damage variables, 01C (1' = I ,II ,III ) are the critical energy release rates associated with fracture modes I, II and III, and (5'! (i=1, 2, 3) are the final displacements corresponding to the complete failure of material respectively. The determination of these parameters is discussed next. 39 Tension 51 . -—> Compressmn N0 damage Damaged Failed (a) Mode I 12.3 A 72,3 .0 (1-d2,3)k2,3 n c f . 62,3 52,3 52,3 GIICJIIC >—> No damage Damaged Failed (b) Mode II and Mode 111 Figure 3.2 Bilinear traction-separation laws 40 3.4.2 Material properties for the cohesive interface The cohesive interface represents a thin (of the order of 10"7 ~ 10"5 mm [27]) resin- rich layer between two adjacent layers, and it can be modeled as an entity of zero (or very small) thickness and large stiffness. The initial interface stiffness does not represent a physically measurable quantity and is introduced as a penalty parameter to simulate a real connection between two adjacent layers before delamination initiation. Ideally, the stiffness of the cohesive element should be infinite so that it does not affect the overall compliance of the model before the damage initiation point. However, a finite value must be used in the finite element context for numerical stability. The initial interface stiffness should be large enough to provide reasonable connections but small enough to avoid numerical problems in finite element analysis. Discussions of different guidelines for selecting interface stiffness can be found in [34]. A penalty parameter of 106 suggested by Zou et al. [33] is used here to multiply the interface strength to calculate the initial interface stiffness. The damage initiation point corresponds to the interface strength, which can be identified as the peak value of the traction-separation relationship. Based on the studies of [28, 34], reducing the maximum interface strength enables obtaining accurate global load-displacement results with coarser mesh. However, the stress concentrations in the bulk material near the crack tip might be less accurate if the used interface strength is too low. The actual interface strength value is used in this work. The damage evolution of cohesive elements defined in terms of traction-separation is controlled by the fracture toughness of the interface (or the critical energy release rate). 41 Fracture toughness values from experiments should be used for the correct prediction of delamination in finite element analysis. 3.4.3 Delamination initiation and propagation criteria The procedure of progressive delamination analysis in ABAQUS/Standard is given in Figure 3.3. The damage initiation and propagation criteria employed in this study will be presented next. 3.4.3.1 Delamination initiation criterion The quadratic stress criteria in terms of tractions of different modes have been shown to predict the delamination initiation successfully in previous studies [35, 28, 33, 31]. Considering that compressive normal tractions do not affect delamination onset, the quadratic failure criterion in two-dimensional is given as [35] II2 +I ’2 I2 :1 (3.1) —c T? 71 IT 2 / where 2'1 is the interlaminar normal stress and 7'2 is the shear stress, Tic (i=1, 2) are the interlaminar normal and shear strengths, and the symbol ( > represents the Macaulay bracket defined by TI 11 > 0 r = 3.2 < l) {0 1'1 < 0 ( ) 42 I Enter load step It I I Calculate stress Is damage initiated? (Eq. 3.1) No Damage —— Is damage propagation criterion met? (Eq. 3.9) Material failed Exit load step It . Crack propagates Figure 3.3 Flowchart of progressive delamination analysis in ABAQUS/Standard Damage evolves 3.4.3.2 Delamination propagation criteria In this section, the energy based delamination propagation criteria are presented. The single-mode delamination propagation criteria are given as follows: Mode I (DCB) £1: (33) GIC Mode II (ENF) 43 .911. = 1 (3,4) GIIC where Gic(i = I ,II ) are the critical energy release rates associated with fracture Modes I and II, which can be determined from experiments. Gi(i = 1,11) refer to the energy release rates, fracture resistance or work done during the damage process by the traction and its conjugate relative displacement in the normal tension (Mode l) and sliding shear (Mode II ) directions, and they are evaluated as G, = 155171ch (3.5) G” = I32 Tzda (3.6) where 5, (i=1,2) are the maximum relative displacements obtained at the point under consideration. The critical energy release rates or the fracture toughness values are equal to the area under the traction-separation displacement curves. That is f G/C =I061 add (3.7) f GIIC =13?- 72615 (3.8) where (fl-f (i=1,2) are the final displacements corresponding to complete the failure of material. The criteria used to predict delamination propagation under mixed-mode conditions are given in terms of the energy release rate and fracture toughness. The criteria available in ABAQUS for the prediction of delamination propagation under mixed-mode 44 conditions are the “power law” by Wu [36] witha =1~ 2, (Eq. (3.9)) and the “BK criterion” (Eq. ( 3.10)) proposed by [37]. a 0: [EL] {ELI =1 (3.9) GIC GIIC GIC +(G/IC -Glc)(-G—) = GC (3.10) T . where 05 =0” in two-dimensional, GT =01 +GS , Gc is the mixed-mode critical energy release rate associated with mixed-mode ratio 2; , 77 is a material parameter T determined by curve fitting. 3.5 Numerical issues of delamination analysis In progressive delamination analyses, which involve nonlinearity due to material softening, the equilibrium equations are solved using a step-by-step incremental solution technique such as the Newton-Raphson method for the complete time range of interest. Because of the linearization assumption used in the Newton-Raphson method, a solution may be subjected to significant errors or may be unstable depending on the time or load step size used. In practice, multiple iterations may be necessary to obtain sufficient accuracy in the solution of the equilibrium equations until the out-of-balance load vector or the displacement increments are sufficiently small. The calculation and factorization of the tangent stiffness matrix can be very expensive for the analysis of large-order systems. One alternative to the Newton-Raphson’s iteration scheme is the Broyden— Fletcher—Goldfarb—Shanno (BFGS) method, which approximates Newton’s method using line search to solve nonlinear problems. The line search technique is activated using the Keyword * Controls in ABAQUS/Standard. 45 -— I: Material softening and stiffness degradation using cohesive elements lead to convergence difficulties in implicit solution procedures like ABAQUS/ Standard . To combat the convergence problems due to material softening, an artificial Duvaut-Lions viscosity model [38] is implemented to cause the tangent stiffness matrix of the softening material to be positive for sufficiently small time increments. Two drawbacks can occur when increasing the viscous parameter. The first one is that using a large viscous parameter may prevent damage localization from happening. Secondly, increasing the viscous parameter will increase the energy dissipated at a material integration point undergoing damage evolution. Consequently, it is important to use as small a viscous parameter as necessary to solve convergence problems [39]. To ensure correct prediction of delamination propagation, the energy dissipated due to viscous regularization must be small compared with the strain energy in the whole model. There are two ways to prevent the delaminated structures from penetrating to each other. One method is to keep the failed cohesive elements. By the constitutive law of the cohesive elements, the stiffness in tension and shear of the failed cohesive elements will be degraded to zero while their stiffness in compression will remain intact. This method is used for the regions where original cohesive elements are being inserted, while contact constraints using keyword “Contact Pair” in ABAQUS are used for the initially cracked regions. 3.6 Model calibration The two-dimensional delamination model is calibrated by comparing the predicted results with the experimental results and numerical results from the 3D model in [31] for the DCB and ENF tests as illustrated in Figure 3.1 (d) and (e) . 46 3.6.1 The specimens and material properties The specimens simulated are 102mm long, 25.4 mm wide, with two 1.56mm-thick arms. The initial delamination length of the specimens is 32.90mm for DCB and 39.30 mm for ENF. Using the notations in Fig. 3.1 (d) and (e), the geometry of the specimens is summarized in Table 3.1, and the material properties used are given in Table 3.2. Table 3.1 Geometry of the specimens [31] (mm DCB I ENF (mm) (mm) 102.00 32.90 | 39.30 1.56 25.40 Table 3.2 Material properties [31] Ply material aroperties AS4/PEEK E11 E22 I E33 I G12 013 I 023 V12 (GPa) (GPa) (GPa) (GPa) (GPa) (GPa) 122.70] 10.10 I 10.10] 5.50 5.50 | 3.70 I025] 0.25 | 0.45 Cohesive layer material properties V13 I V23 Stiffness (Calculated) Strength Fracture toughness (MPa) (MPa) ( J /m2) 1‘1 I k2 I k3 if I r; I 736 GIC I GIIC I Gmc ii 80X106 1100x106 IN/A 80 I 100 IN/A 969 I1719 IN/A 3.6.2 Selection of elements and mesh size prediction The incompatible mode plain strain elements CPE41 from ABAQUS [l] are selected for the plies and the two-dimensional cohesive elements COH2D4 are selected for modeling the zero-thickness interface layers. The incompatible mode elements CPE41 have incompatible deformation modes added in addition to the standard displacement degrees of freedom to eliminate the parasitic shear stresses that cause the response of the 47 regular first-order displacement elements to be too stiff in bending. Furthermore, they use full integration, and thus they have no hourglass modes. Different guidelines and formula for the estimation of the cohesive zone length and the maximum mesh size in cohesive zone in the direction of crack propagation have been discussed in the literature [34, 40]. Turon et al. [34] proposed a formula to predict the cohesive zone lcz length as follows GIC (If)2 where E33 is the Young Modulus of the material associated with Mode 1, GIC is the (3.11) lcz = M533 Mode I fracture toughness, if is the maximum interfacial strength, and M is a parameter that depends on each cohesive zone model. Following [41] (M =1) and applying the material properties from Table 3.2 in Eq. (3.11), the cohesive zone length is calculated to be 1.53mm. According to the study in [42], a minimum of three elements should be included in the cohesive length to capture the smoothing continuum of randomly inhomogeneous material so that each of the three regions of the cohesive zone i.e. elastic, damage evolution, and fully damaged, can be represented by one element. By allowing three elements in the cohesive zone along the direction of crack propagation [34,42], the element size in the direction of crack propagation is given as _ 1.53mm 16 = 0.51mm The mesh converge study with different mesh sizes 1m, 0.5m, 0.25mm, and 0.125mm using a DCB test is given in Figure 3.4; and this validates the element size estimation from the calculations above. The study shows that a mesh size of 0.5mm based 48 on the given material properties gives converged results, and it is used for all the following two-dimensional case studies of progressive delamination in this chapter. Figure 3.5 shows the stress distribution at the delamination front from a DCB test. 200. - . ~ . , 150- , A 2 v § 100 - . o 15 so» MESH_0.125mrn . MESH_0.25mm — MESH_0.5mm -— MESH_1.0mm 0. . r . 1 r r A 1 . r 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Displacement(mm) Figure 3.4 Cohesive element mesh converge study using a DCB test 49 S, 511 (Avg: 75%) +4.574e-1-08 +3.811e-1-08 7 +3.049e+08 +2.287e+08 _ +1.525e+08 +7.625e-1-07 +3.181e+04 -7.619e+07 -1.524e+08 -Z.286e+08 -3.049e+08 -3.8118+08 -4.573e+08 (a) Tensile stress in the plies S, 522 (Avg: 75%) +5.9066+07 +5.2806-I-07 +4.654e+07 +4.027e-1-07 +3 .401e+07 +2775e+07 ” +2.149e+07 +1 .5226-1-07 +8.960e+06 +2.697e+06 '3. 566e+06 -9.829e+06 -1.609e+07 (b) Nominal inter-laminar normal stress in the cohesive zone Figure 3.5 Stress distribution at the delamination front in a DCB test 3.6.3 Comparison of results Using the elements and mesh size selected and single mode delamination propagation criteria, DCB and ENF tests in [31] are simulated and the comparisons of the results are 50 given in Figure 3.6 and in Table 3.3. It can be concluded that good agreements have been obtained between the current two-dimensional CZM and the 3D numerical model and experimental results in [31]. Table 3.3 Summary of result comparison of model calibration Tests Max. Load (N) Displacement at Max. Load Experiment Numerical Experiment Numerical [31] 30131] ICurrentZD [31] 30131] I Current2D DCB 147.10 135.30 | 138.50 4.80 4.70 | 4.50 ENF 734.00 696.00 | 720.50 3.90 3.90 | 3.80 3.7 Selection of mixed-mode delamination propagation criterion To determine the appropriate delamination propagation criterion, the mixed-mode criteria presented earlier are studied and compared with the single mode criterion under Mode I loading condition using an antisymmetric cross-ply [0/90] double cantilever beam. The results given in Table 3.4 demonstrate the existence of mix-mode effect in [0/90] laminate under single mode (Mode 1) loading condition. The results show that the forces required to peel the cross-ply double cantilever beam using different delamination propagation criteria are in the following order: Single mode < BK law < Power law (a = 1) \ (a) [0/0] (b)[0/90] (c)[90/90] Figure 3.9 Post-test shapes of 2 ply laminates in DCB 56 3.8.2 DCB test of four-ply laminates For a comparison of all the laminates studied, the load-displacement curves obtained from the DCB test are given in Figure 3.10. The unsymmetric laminate [0/90/0/90]t gives slightly higher delamination resistance than the symmetric laminate because Mode 11 is involved in the unsymmetric laminate, which has coupling between extension and bending. Since the load-displacement curve for [0/90]s is almost indistinguishable from that of [90/0]s, only result of [0/90]s is plotted in Figure 3.10. Similar to the conclusions from the two-ply laminate study, the study using the four-ply laminates shows that the stiffness is a dominant factor in the Mode 1 delamination resistance. 3.8.3 Comparison of plain weaves with 4-ply laminates in DCB In this section, composites made of symmetric plain weaves of [0-90/90-0]s and [90-0/0- 90]s, and antisymmetric plain weave [0-90/90-0/90-0/0-90]t as shown Figure 3.11 are simulated in DCB tests. The results are compared with their cross-ply laminated counterparts as shown in Figure 3.12. As can be seen in the figures, the load- displacement curves corresponding to the laminated regions match very well between the plain weave composites and the associated cross-ply laminates. However, the load drops to a plateau segment through the undulation region in the plain weave because the average stiffness is lower in the undulation region. Similar to the conclusion of the 4-ply laminate study, the antisymmetric laminate of plain weave gives slightly higher delamination resistance than the symmetric laminates because Mode II is involved in the nonsymmetric laminates due to the presence of the coupling between extension and bending. 57 Forum) 150. -— DCB [0/0] ’ — ocetolso] --------- oce [woo/0190]: ' DCB [O/90]s 100. r ---- DCB[90/90] I/{AM'M-E‘Kh '63). Why,“ 2 50' ,_ "’ MM “Mm-u... _ I,“ ........... «mmh'mcm and.» ,c' ........................ ..h ............... ’1’ ------------------- :v ............... .3"; ......... / ...... o C: ...... 0. 5. ‘ 10. 15. ‘ 20. ‘ 25. o. Displacemenqmm) Figure 3.10 DCB load-displacement curves of laminates 0 0 90 90 (a) [0-90/90-0]s 0 0 (b) [90-0/0-90]s R90 0 (c) [0-90/90-0/90-0/0-90]t Figure 3.11 4-ply plain weave models 58 150. . , . , - , . , . . . -- PW_[0-90/90-0/0-90/90—0]t — PW_[0-90/90-0]s -- PW_[90-0/0-90]s ----- lam_[0/90/0/90]t ----- lam_[0/90]s ----- lam_[90/0]s 100. - Force(N) 50 . '- Mth-m d 0. . 1 . 1 . 1 r 1 . r . o. 5. 10. 7 15. 20. 25. 0. D1splacement(mm) Figure 3.12 Comparison of DCB load-displacement curves of plain weaves and laminates 3.8.4 Undulation effects on delamination propagation in plain woven To study the effects of the undulations in woven composites on the delamination propagation, three models of the symmetric plain weave [90-0/0-90]s as illustrated in Figure 3.13 are analyzed in this section. In Model A (Figure 3.13 (a)), a single “horizontal” cohesive layer is inserted into the middle plane of the laminate to simulate the matrix interface at the middle plane and restrict the delamination to propagate along this interface; In Model B (Figure 3.13 (b)), a single “skew” cohesive layer is inserted above the upper surface of the 0 ply of the lower arm through the undulation region to confine the propagation of the delamination along the matrix interface there; In Model C (Figure 3.13 (c)), both the “horizontal” and “skew” cohesive interfaces in Mode A and Model B are inserted to allow the delamination to propagate along both paths simultaneously. 59 11111 I W“ Wotan" I III "IIIIIIIOIII." nursunumsawurmw W :1:in meat 1. ...-.1 1.1. 211 r I r u s urrlulnnuuuumuuur "W m- rerun-0.1mm .. w :4 mm. 1ch tin miliuir miimmt;l;nu 90 90 0 muwmllIllllllflllllllllllllllllllmlillilllllllllllllIlIIIIIIIIlIIIIIIIIIIII IIIIIIIIIIIIIIIIIIIIIIIIIIIIlllllltmlllltluulmttluullllm1:11 111111. 18111011111111111111111111111a. 11111 11111181111111“ 111111111 I lIIIIIIIIII"IIIIIIIIIIIIIIIIIIIIummmummlMlIl-m ”11111111111 11111111111. mu 111 1:1 111 111111- 0131311111 11. mm 1111 I' "‘hllmm1mmumnmmuimmumummmumnm1111mm (b) Model B with “skew” cohesive layer (c) Model C with both “horizontal” and “skew” cohesive layers Figure 3.13 Three models of plain weave [90-0/0-90]s 3.8.4.1 DCB test results and discussions The load-displacement curves of the three models from the DCB tests are given in Figure 3.14. In the legend, “Horizontal” denotes Model A, “Skew” denotes Model B, and “Hor_Skew” denotes Model C. The three load-displacement curves are almost identical before the undulation region. Once the delamination reaches the undulation region, the three curves start to digress from each other. Model A with the horizontal cohesive layer stays on the top, which means it consumes more energy for the delamination to propagate along the horizontal interface. Model B with the skew cohesive layer stays on the bottom, which means it consumes the least energy for delamination to propagate along the skew interface. The load-displacement curve of Model C with two cohesive layers shows some fluctuation when the delamination enters the undulation region, and then it drops to 60 match the curve in Model B. The fluctuation phenomenon or the so-called bifurcation effect in Model C happens because of the competition of the two cohesive elements at the intersection between the cohesive layers. According to [44], the bifurcation effect is a numerical problem inherent in a softening material model and it happens when multiple cohesive elements experience damage simultaneously. Eventually, the delamination in Model C propagates along the skew interface as in Model B because it consumes less energy along this path. This is shown by the load-displacement curve overlap between Model B and Model C. The DCB simulation results may be explained by comparing the stiffness of the materials on the two sides of the interfaces. The material distribution is symmetric about the horizontal cohesive interface, and thus the materials have the same stiffness on the two sides of the horizontal interface. However, the material on one side of the skew layer is different from the other side; and thus, a stiffness mismatch exists between the two sides of the skew interface. This stiffness mismatch may lead to the easy delamination propagation along the skew interface. 3.8.4.2 ENF test results and discussions In this section, the three models are simulated in ENF tests. The resulted load- displacement curves of the three models are given in Figure 3.15. Contrary to the delamination propagation in the DCB tests, the ENF results show that it consumes less energy for the delamination due to Mode II to propagate along the horizontal cohesive layer than along the skew cohesive interface. This may be explained by the elementary beam theory, according to which, the peak shear stress occurs in the 61 100. - 80$ — Horizontal — Skew 1 - ------- Horizontal and Skew g 50.1 1 1 2 11’. 40.. l 20.~ 0. . T o. 5_ , 15. 2111. Displacement (mm) Figure 3.14 Comparison of DCB Load-displacement curves using the three models Force (N) (l and 2 denote the undulation locations in Figure 3.13 (a)) 500. * 400.~ 1 i 300 * 200 1 — Horizontal 100.- _ Skew ------- Horizontal and Skew 0. ‘ I L l (11. 4 8° ' Displacement (mm) 12. Figure 3.15 Comparison of ENF Load-displacement curves using the three models (“1” denotes the undulation location in Figure 3.13 (a)) 62 middle plane of the beam. Again, the bifurcation effect is observed in Model C with two cohesive layers as shown in the load-displacement curve in Figure 3.15. 3.8.4.3 Energy dissipation due to viscous regularization To ensure correct prediction of delamination propagation, the energy dissipated due to numerical viscous regularization must be small compared with the strain energy in the model. As shown in Figure 3.16, the energy dissipation due to numerical viscous regularization is zero in the DCB tests for both Model A and Model B. In the ENF tests, the energy dissipated due to numerical viscous regularization is less than 5% of the strain energy for both Model A and Model B as illustrated in Figure 3.17. This shows that the numerical viscous regularization and stabilization schemes used in the analyses do not affect the correct prediction of the final results. The presence of the dissipation energy in ENF tests may be associated with the stabilization scheme used to prevent penetration between upper and lower beams in the initial crack area. 63 0.80 0.60- Energy(J) p 8 0'20" ...... Dissipation Enargy(Horizontal) _ — Strain EnergyO-lorizontal) 0.00 1 - 1 - C1. 5. 10. 15. IO. Displacemenflmm) (a) Model A 0.80 a , - 1 0.60— S >. 2' 0.40 o c ill 0'20" ''''''' Dissipation Enorgy($kaw) -- Strain Enorgy(Skow) 0.00 1 - 1 01. 5. 10. 15. O. Displacemenflmm) (b) Model B Figure 3.16 Energy dissipation in DCB tests 64 2.5 2.0- 0.5 2.5 J ........ Dissipation Energy (Horizontal) -— Strain Energy(Horizontal) J A. ' Dlsplacemenflmm) ' (a) Model A 2.0- 0.5 ------- Dissipation Enargv($k¢w) — Strain Energy(Sk°W) L A J A J 4. B. 12. Displacemenqmm) (b) Model B Figure 3.17 Energy dissipation in ENF tests 65 3.8.4.4 ENF test results and discussions To further study the undulation effects on delamination propagation in woven composites, the whole undulation curves in plain woven composites are simulated with cohesive layers added into different interfaces in the models. The initial crack and the overall geometries of the specimen are the same as given in Table 3.1; and the dimensions of the undulation regions are illustrated in Figure 3.18. Model D, E, and F (See Figure 3.19) with cohesive layers simulating the upper, middle, and lower delamination paths in the undulation regions are studied by simulating ENF tests. Figure 3.20 (b) shows the resulted load-displacement curves. In Figure 3.20 (a), the numbers denote the location where the undulations occur in the plain weaves. As can be seen from Figure 3.20 (b), delamination propagation through the undulated paths in Model D and F consumes more energy than that through the horizontal layer in Model B. This conclusion is the same as that from the study of models with half undulation path in plain woven composites. The load-displacement curve of Model D with upper cohesive layer shows smooth transition through undulation “2” and “3” compared to the other two models. This may be due to the difference in the mesh used: 0.25mm is used for the cohesive layer for Model D due to convergence issues, while 0.5mm is used for the other two models. V A F V 1 93.4, .u - . .. ‘33:}; - ...5' 05.1.1 -,i ' tz‘|"-':f1'i’.‘-"';.;::aivstall!Qiiiifi'iéiis Ei’E"‘i§’.di".ii"'." 1‘ ' i ‘ ' ' - -- . 'a “hunk jam" mag: "'.= ._ z 3 12 - " :1 ‘:1 . '3. 51.1..“ 91,- : '31:: .:'-"1 ' 1 ‘ . ‘ ‘ . ~13. 1. 1‘ .. . . a 'a‘ - ‘1 ‘ ' A IV 7" Figure 3.18 Dimensions of the undulation region (Dimensions in mm) 66 (a) Model D with upper cohesive layer (b) Model B with horizontal cohesive layer Mtétunahhnfinhv-uwu - - = ~' mmumumumm IIIlllllllllIIIIIIIIIIIIIIIIIIIIIIIIIIIIIltllulllllll"mm “iglltlllillillilllflllllmlum IIIIIIIIIIIIIIIIIWIIIIMIIIMIIIIIIIII ummmmnw III! nun-urn- ‘fG.-'€* mull] umnmnmmmwmpnl (a... c. I (0) Model F with lower layer Figure 3.19 Models with through cohesive layers '31! >axdbli>53t Li. ‘ v. ‘ -: 1 . . ' ' . .’ '.“ ‘3..th1‘41«fri{tkhfi:§l‘f{£*ftrxixi3&i 600. f . . 4 500.~ 3 . . 2 g..' 400. ~ '1 . f l \ 1 ,. . §3oo¢ - o \ l l 200. " —- Upper " '— Middle * 100. Lower 0. 5. 10. 15. 20. Displacementtmm) (b) Load—displacement curves of Model D, E and F with numbers denoting the undulation locations Figure 3.20 Load-displacement curves of Model D, E and F 67 For comparison purpose, the displacement curves of Model A (Figure 3.13 (a))with horizontal cohesive layer in plain weave model of one-undulation and Model E (Figure 3.19 (b)) with horizontal cohesive layer in plain weave model of two undulations are depicted in Figure 3.21. As can be seen, the peak loads are the same for both models. The difference in the displacements at which the peak loads occur is due to the stiffness difference in the two models: model B has lower stiffness than Model A because of the one additional undulation involved. 600. 500. ~ .... C ........ -- Middle (one undulation) 100' . """ Middle (two undulations) - o. ' 5'. ' '16. ' 1'5. ' 2’0. Displacemenflmm) Figure 3.21 Load-displacement curves of Model A (one undulation) and E (two undulations) 3.9 Conclusions of two-dimensional delamination study 0 In laminated composites, the average stiffness of each arm of a DCB is an important factor for delamination initiation and propagation. Large stiffness tends 68 to increase the maximum load required for delamination initiation, but it also promotes the delamination propagation rate once delamination propagation starts. 0 Compared to cross-ply laminated counterparts, plain weave composites have a plateau in their DCB load-displacement curves associated with the undulation region because of the stiffness variation from the laminated regions. 0 Nonsymmetric laminates give slightly higher delamination resistance than their corresponding symmetric laminated counterparts because Mode II is involved in the nonsymmetric laminates due to the presence of the coupling between extension and bending. 0 The delamination propagation direction through the undulation regions in plain woven composites depends on the loading conditions: 0 Delamination due to Mode I in DCB test propagates along the skew interface through the undulation region. 0 Delamination due to Mode II in ENF test propagates along the skew interface through the undulation region. Based on this result, delamination resistance may be improved by eliminating the horizontal pure matrix interface in laminated composites using interlocking through layers as in the Q3D weave design. 0 Cohesive layers should be used selectively due to the bifurcation effect. 0 The two-dimensional study is computationally efficient, but the information that can be obtained is very limited compared to 3D analysis. 3D progressive damage analysis with intra-laminar damage modes will be considered in the next chapter. 69 10. 11. 12. REFERENCES ABAQUS/Standard User’s Manual, online ver. 6.8 (2008), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. C. T. Herakovich, Mechanics of fibrous composites, John Wiley & Sons, Inc. , 1998 A.J. Brunner, B.R.K. Blackman, and P. Davies, A status report on delamination resistance testing of polymer-matrix composites, Engineering Fracture Mechanics 75 (2008) 2779 - 2794 AC. Orifici, I. Herszberg, and RS Thomson, Review of methodologies for composite material modelling incorporating failure, Composite Structures, Volume 86, Issues 1-3, November 2008, Pages 194-210 T.E. Tay, Characterization and analysis of delamination fracture in composites: An overview of developments from 1990 to 2001, Appl Mech Rev vol 56, no 1, January 2003. D. Liu, Impact-Induced Delamination—A View of Bending Stiffness Mismatching, Journal of Composite Materials, Vol. 22, No. 7(1988), 674-692 P. Naghipour, M. Bartsch, L. Chemova, J. Hausmann, and H. Voggenreiter, Effect of fiber angle orientation and stacking sequence on mixed mode fracture toughness of carbon fiber reinforced plastics: Numerical and experimental investigations , Materials Science and Engineering A 527 (2010) 509 - 517 J. Andersons and M. Konig, Dependence of fracture toughness of composite laminates on interface ply orientations and delamination growth direction, CompOs Sci Technol 64 (2004). PP. 2139—2152. N .S. Choi, A.J. Kinloch, and JG Williams, Delamination fracture of multidirectional carbon—fiber/epoxy composites under mode 1, mode 11 and mixed-mode I/II loading, J. Compos. Mater. 33(1999), 73-100. J. Tao and CT. Sun, Influence of ply orientation on delamination in composite laminates, J. Compos. Mater. 32 (1998), 1933—1947. A.B. de Morais, M.F. de Moura, J .P.M. Goncalves, and PP. Camanho, Analysis of crack propagation in double cantilever beam tests of multidirectional laminates, Mechanics of Materials 35 (2003) 641—652 A.P. Mouritz, C. Baini and I. Herszberg, Mode I interlaminar fracture toughness properties of advanced textile fibreglass composites, Compos Part A: Appl Sci A20 (7) (1999), Pp. 859-870 70 13. Y. Tanzawa, N. Watanabe, and T. Ishikawa, Interlaminar fracture toughness of 3-D orthogonal interlocked fabric composites, Compos Sci Technol 59 (8) ( 1999), pp. 1261—1270. 14. J. Chen, E. Ravey, S. Hallett, M. Wisnomc, and M. Grassi, Prediction of delamination in braided composite T-piece specimens, Composites Science and Technology 69 (2009) 2363—2367 15. A. E. Bogdanovich, “Multi-scale modeling, stress and failure analyses of 3-D woven composites” J Mater Sci (2006) 41:6547—6590 16. X. Tang, J. D. Whitcomb, A. D. Kelkar, and J. Tate, Progressive failure analysis of 2 x 2 braided composites exhibiting multiscale heterogeneity, Composites Science and Technology 66 (2006) 2580 - 2590 17. S.V. Lomov, DS Ivanov, I. Verpoest, M. Zako, T. Kurashiki, and H Nakai, Meso-FE modelling of textile composites: Road map, data flow and algorithms. Compos Sci Technol 2007; 67 :1870 - 91. 18. L. Gorbatikh, D. Ivanov, S. Lomov, and I. Verpoest, On modelling of damage evolution in textile composites on meso-level via pr0perty degradation approach, Composites: Part A 38 (2007) 2433-2442 19. R. Krueger, Virtual crack closure technique: history, approach and applications. Appl Mech Rev 2004;57:109—43. 20. 1.8. Raju, Calculation of strain-energy release rates with higher order and singular finite elements, Eng. Fract. Mech. (1987), 28(3), 251-274. 21. E.F. Rybicki and M.F. Kanninen, A finite element calculation of stress intensity factors by a modified crack closure integral, Eng. Fract. Mech. (1977), 9, 931—938. 22. C. Balzani and W. Wagner, An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates, Engineering Fracture Mechanics, Vol. 75, Issue 9, June 2008 2597-2615 23. W. Wagner, F. Gruttmann, and W. Sprenger, A finite element formulation for the simulation of propagating deleminations in layered composite structures, Int J Numer Methods Engng 51 (2001). Pp. 1337-1359 24.2. Petrossian and M. R. Wisnom, Prediction of delamination initiation and growth from discontinuous plies using interface elements, Composires Part A 29A (1998) 503-515 25. A. Pantano and R. C. Averill, A mesh-independent interface technology for simulation of mixed-mode delamination growth, International Journal of Solids and Structures Vol. 41, Issue 14, July 2004 3809-3831 71 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. RP. Camanho, and CG. Davila, Mixed—Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials, NASA/1‘ M-2002-21 1737 Z. Zou, S.R. Reid, S. Li, and RD. Soden, Application of a delamination model to laminated composite structures, Composite Structures 56 (2002) 375—3 89 G. Alfano and M. A. Crisfield, Finite element interface models for the delaminationanalysis of laminated composites: mechanical and computational issues, Int. J. Numer. Meth. Engng 2001; 50:1701-1736 G.I. Barenblatt, Mathematical Theory of Equilibrium Cracks in Brittle Failure, Advances in Applied Mechanics 7 (1962). A. Needleman, A Continuum Model for Void Nucleation by Inclusion Debonding, Journal of Applied Mechanics 54 (1987) 525—53 1. PP. Camanho and CG. Davila, Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials, NASA/TM-2002-211737 O. Allix and P. Ladeveze, Interlaminar interface modelling for the prediction of delamination. Compos. Struct. 22 (1992) 235-242. Z. Zou, S.R. Reid, and S. Li, A continuum damage model for delaminations in laminated composites, Journal of the Mechanics and Physics of Solids 51 (2003) 333 — 356 , A. Turon, C.G. Davila, P.P. Camanho, and J. Costa, An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models, Engineering Fracture Mechanics 74 (2007) 1665—1682. W. Cui, M. R. Wisnom, and M. Jones. ”A Comparison of Failure Criteria to Predict Delamination of Unidirectional Glass/Epoxy Specimens Waisted Through the Thickness.” Composites 23(3), 1992, 158-66 E.M. Wu and Jr C. Reuter, Crack extension in fiberglass reinforced plastics, Report no. 275, T&AM, University of Illinois, 1965. ML. Benzeggagh and M. Kenane, Measurement of Mixed-Mode DelaminationFracture Toughness of Unidirectional Glass/Epoxy Composites With Mixed-Mode Bending Apparatus, Composites Science and Technology 56 (1996) 439—449. G. Duvaut, Lions, J .L., Inequalities in Mechanics and Physics, Springer-Verlag Berlin Heidelberg New York, 1976 P. Maimi, PP Camanho, J Mayugo, and CG Davila, A continuum damage model for composite laminates - part II: computational implementation and validation.Mech Mater 39 (2007) 909—19 72 40. P. W. Harper and S. R. Hallett, Cohesive zone length in numerical simulations of composite delamination, Engineering Fracture Mechanics 75 (2008) 4774—4792 41. A. Hillerborg, M. Modéer, and RE. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cement Concr Res 6 (1976), pp. 773-782. 42. Z.P. Bazant and BH. Oh, Crack band theory for fracture of concrete, Materials and structures 1983 (16) 155—177 43. ,J. R. Reeder ”An Evaluation of Mixed-Mode Delamination Failure Criteria.” NASA TM 104210, 1992 44. K. Y. Volokh, Nonlinear Elasticity for Modeling Fracture of Isotropic Brittle Solids, J. Appl. Mech. ,January 2004 -- Volume 71, Issue 1, 141 73 Chapter 4 Three-dimensional static progressive damage study of laminates 4.1 Abstract This chapter is aimed to show the necessity of in-ply matrix cracking for the prediction of delamination in laminated plates subjected to transverse static loading. First, a three-dimensional continuum damage model is developed for the prediction of the initiation and propagation of intra-laminar damage mechanism. This continuum damage model is then combined with the cohesive zone model in ABAQUS/Standard [1] to perform three-dimensional progressive damage analysis of skew-angled laminated plates subjected to transverse static loadings. Finally, the numerical predictions are compared with the results of the delamination patterns and fiber-bridging damage modes from the original experimental works in [2,3]. Good correlations are achieved between numerical simulation and experiments for the damage modes and delamination patterns. 4.2 Literature review and the modeling methodology 4.2.1 Continuum damage mechanics models for composites In Continuum Damage Mechanics (CDM) models e.g. [4], material damage, the initiation and growth of micro-cracks in a continuous medium was modeled with internal damage variables first introduced by L.M. Kachanov [5]. The damage variables can be interpreted as the reduction of effective loading area to characterize the degradation of the material stiffness caused by damage such as cracks and voids. Effective stresses were introduced by Y.N. Rabotnov [6] to represent the stresses acting on the net (undamaged) loading area. The effective stresses may be expressed in terms of the damage variables and the actual stresses by applying equilibrium equivalence in one-dimension (1D) model 74 and energy equivalence in three-dimensional model as explained later. The effective stresses will be used in the failure criteria to define the loading functions later. Crack-bl A0 (a) Damaged material with area A0 under actual tress 0' 6 Ad=(1—d)Ao 6 (b) Damage model with net loading carrying area Ad under effective stress 0'" Figure 4.1 One-dimensional continuum damage mechanics model To introduce the concepts of damage variables and effective stresses in CDM, a one- dimensional CDM model is given in Figure 4.1. In Figure 4.1(a), damage is created in the form of a crack by an actual stress 0 in a specimen with cross section area/10. The continuum damage model for this one-dimensional example is given in Figure 4.1(b), where Ad denotes net (undamaged) loading area in the fictitious undamaged configuration, and 6' is the effective stress applied on this fictitious area. To characterize the degradation of material stiffness due to damage, damage variable (1 (OS (1 $1) is introduced and defined as follows. d :10? (4.1) Eq. (4.1) can be rewritten as follows: 75 Ad = (1— 60/10 (4.2) which reveals that the net loading carrying area Ad that characterizes the stiffness of the material will decrease with damage growth. A zero value of damage variable indicates undamaged material, and a value of 1 represents fully damaged material. Applying equilibrium equivalence to the actual configuration Figure 4.1(a) and fictitious undamaged configuration Figure 4.1(b), we have 6Ad = 0A0 (4.3) Inserting Eq. (4.2) into Eq. (4.3), the effective stress can then be defined in terms of the damage variable (I and the actual stress 0 as follows. 6' = 3— (4.4) l — d The effective stress in the fictitious undamaged configuration is imaginary and can not be measured directly, while the actual stress a is a physically measurable quantity. As can be seen in Eq. (4.4), the effective stress increases as damage accumulates in the specimen. The so-called damage operator 1 is a magnifying factor that relates the actual stress to the effective stress. In the one-dimensional case, the damage variable and effective stress are scalars. As shown later, they are tensors in three-dimensional CDM models. Assuming orthogonal damage modes and using the energy equivalence principle, which assumes that the elastic energy of the damaged material is in the same form as that of the undamaged material except that the actual stresses are replaced by the effective stresses; Cordebois and Sidoroff [7] outlined a three-dimensional damage model to take 76 into account the material anisotropy. Matzenmiller [8] presented a detailed formulation of constitutive model for fiber-reinforced composites with elastic-brittle behavior. Based on the anisotropic damage mechanics models developed by Cordebois and Matzenmiller, Zako [9] performed a three-dimensional progressive damage analysis on a unit cell of plain woven fabric composites subjected to on-axis uniaxial loading using stress-based failure criteria. Starting from the definition of complementary free energy density, Maimi et al. [10- 12] developed a continuum damage model for the prediction of the onset and evolution of intra-laminar failure in composite laminates under plane stress conditions. Using equivalent displacement and stress concepts, Lapczyk et al. [13] presented an anisotropic damage model for predicting failure and post-failure behavior in composite laminates under plane stress. To alleviate mesh dependency due to strain softening, both works used the crack band model to integrate a characteristic length into the constitutive models to calculate the strain corresponding to the material failure. 4.2.2 Progressive damage analysis of laminated plates subjected to transverse loading Using CZM for both the intralaminar and interlaminar damages, the peanut shape delamination in cross-ply laminates was simulated in [14,15]. But this approach requires the knowledge of the locations of the intra-laminar damage in advance. Using CZM for delamination and stress-based criteria for interlaminar damage, the peanut-shaped delamination in cross-ply laminates subjected to transverse loads is reproduced in [16, 17]. In their intralaminar damage models, once the stress-based criterion for an in-ply damage mode is met, the associated elastic modudi or stresses at 77 corresponding material points are reduced or zeroed accordingly. This method needs pre- runs to determine the parameters describing the damage-resisting capability of the structures [18]. Most recently, by extending the continuum damage model introduced by Maimi et al. into a three-dimensional model for intralaminar failure prediction and using the cohesive zone model for interlaminar failure prediction, Lopes et al. [19] performed a progressive damage analysis of dispersed stacking sequence laminates subjected to transient low velocity impact using the explicit finite element method in ABAQUS. Pilchak et al. [3] investigated experimentally the fiber-bridging damage modes in small-angle laminates and proposed a technique to combat delamination using small- angle laminates. To the author’s best knowledge, simulation of the fiber-bridging damage modes in small-angle laminates subjected to transverse loads has not been reported in the literature. 4.3 Three-dimensional progressive damage modeling methodology in this work The CZM approach is efficient for the simulation of delamination, a fracture at the interface in composites. It is less efficient for modeling intralaminar damage, the fracture in the bulk materials results in the simultaneous use of two material models for the same physical material. It is desired to incorporate a material failure law directly in the stress- strain constitutive description of the bulk material. This is achieved by modeling the material failure via elastic strain softening model with stiffness degradation to be introduced in the next section. The overall three-dimensional progressive damage modeling methodology adopted here is shown in Figure 4.2. 78 Progressive Analysis Model for Fiber-Reinforced Composites Y i Intra-laminar Damage (Fiber damage, matrix Inter-laminar Damage (Delamination) cracking/Crushing) Continuum Damage Cohesive Zone Model Mechanics Model (Camanho 2002 ) i l i Damage Initiation Damage Propagation (Stressed Based Failure (Fracture Mechanics Criteria) Based Failure Criteria) l Mesh Sensitivity Alleviation Using Crack Band Model (Bazant1983) Figure 4.2 Three-dimensional progressive damage modeling methodology 4.4 Development of the continuum damage mechanics model for intra-laminar damage In this section, a continuum damage model for intra-laminar damage is developed. The basic concepts, constitutive models, and numerical techniques to be introduced about this model are given in Figure 4.3. 79 Continuum Damage Mechanics Model i v i Basic Concepts: Constitutive Model: Numerical Techniques: - Damage modes - Energy dissipation - Crack band model - Damage variable - Consistency conditions - Viscous regularization - Effective stress - Damage evolution laws - Damage operator - Damage threshold - Loading function - Activation function Figure 4.3 Continuum damage model 4.4.1 Formulation of one-dimensional continuum damage model This section is aimed to introduce the basic concepts and formulation of a continuum damage model using a uniaxial example. 4.4.1.1 The elastic softening model Fiber-reinforced polymer composites behave like quasi-brittle materials, i.e., such materials fail after reaching peak strengths without significant plastic deformation [8]. Consequently, fiber-reinforced polymer composites can be modeled as elastic-brittle materials with neglected plastic deformation. In the elastic softening model with stiffness degradation [20], strain-softening is due solely to the degradation of material stiffness without other inelastic behavior. After the material is damaged, unloading leads to the complete closure of the cracks. This means that the material will unload along a straight line pointed toward the origin in the stress-strain curve of the material and reload at the degraded stiffness along the line with a slope of (l-d)E as illustrated in Figure 4.4. 80 In Figure 4.4, 0'6 is the strength of the material and EC is the strain associated with 0'6 at which the material damage (strain softening) initiates and can be calculated as 0C EC = —— (4.5) E where E is the Young’s Modulus of the material. The damage variable (1 represents the degradation of material stiffness. 8f is the strain associated with material failure. Then the questions to be answered in Figure 4.4 are: How is 8f computed? How is the damage variable (1 calculated? The first question will be answered by introducing the concept of mesh dependency and crack band model, and the second question will be answered later. 0" Figure 4.4 Elastic softening model with stiffness degradation 4.4.1.2 Mesh dependency and crack band model In continuum mechanics, the constitutive model is expressed in terms of stress-strain relations. Once the material is damaged, the material exhibits strain-softening behavior leading to strain localization. The softening stress-strain constitutive equations may result 81 in physically inadmissible responses. That is, the energy dissipated is proportional to the volume of the failed element. As a result, the damage is localized in a region of zero volume, which is the so-called softening crack band [20] as shown in Figure 4.5. This leads to a fracture without energy dissipation with mesh refinement. This is known as mesh dependency of the finite element results in that the energy dissipated decreases upon mesh refinement. A good example of mesh dependency can be found in [13]. To correct this behavior, the strain associated with material failure is computed using the crank band model proposed by [21], in which fracture is modeled as a smeared crack band consisting of parallel densely distributed micro-cracks as illustrated in Figure 4.5 0 4 Softenin- crack band 0 Figure 4.5 Stress-strain curve of the crank band model In Figure 4.5, we is the crack band width or the so-called characteristic length in the unloaded state, a strain softening localization limiter to overcome mesh-dependency, A is the elongation of the crack band in loaded state, — 1S the inelastic fracturing straln, WC and% is the elastic strain of the bulk material. Using the crack band model, 8f is calculated by keeping the energy dissipated over the fracture surface constant as follows: 82 20‘” C awe sf = (4.6) where CF is the fracture toughness. 4.4.1.3 Damage variable calculation Once the material is damaged, the damage evolution is controlled by the damage variable defined as follows: f max_ c (1:8 (8 8) (4.7) gmaxwf _8c) x where am is the maximum strain obtained after damage initiation and before full damage, and it is used to track the damage in the strain domain and can be calculated as follows. emax =min{ rnax {86,8},€f} (4.8) i =0” max Eq. (4.8) indicates that the maximum strain emax must increase monotonically after damage initiation associated with strain EC and before full damage associated with failure strain 8f . It specifies the condition that material can not be recovered after damage, which ensures thermodynamics irreversibility. To scale the strains by the critical strain 8C associated with material damage initiation, the elastic domain threshold is defined as follows. r = — (4.9) As can be seen in Eq. (4.9), the elastic domain threshold r is a dimensionless scalar. Using Eq (4.9), the damage evolution can be mapped from strain domain as shown in Figure 4.4 to the domain of elastic domain threshold as illustrated in Figure 4.6. 83 ‘Y 0 Figure 4.6 Damage evolution in the domain of elastic domain threshold The damage variable can then be expressed in the elastic domain threshold as f max_ - r (r 1) (4.10) — rmax (rf —1) . . . . . . . max . where rf IS the elastic domam threshold assoc1ated w1th material failure, and r 1S the elastic domain threshold associated with the maximum elastic strains that can be attained before further damage occurs, i.e. rf =—— (4.11) and rm" = 8 (4.12) Using Eqs. (4.8), (4.11), and (4.12), we have rmaxzmin{ max {1,r},rf} (4.13) t=0,tmax Eq.(4. 13) specifies the thermodynamics irreversibility condition as Eq. (4.8) does. 84 4.4.2 Formulation of three-dimensional continuum damage model After the introduction of the basic concepts and formulation of the one-dimensional continuum damage model, the derivation of the three-dimensional orthotropic continuum damage model is presented in this section. 4.4.2.1 Damage modes in fiber-reinforced polymer composites Dominant damage modes in a fiber-reinforced polymer composite subjected to transverse loading include intralaminar damage e.g. fiber breakage, matrix cracking, matrix crushing, and interlaminar damage (delamination). The fibers are in forms of bundles and arranged unidirectionally within lamina of laminated fiber-reinforced polymer composites; and thus, they may be treated as orthotropic materials with longitudinal fiber damage mode and transverse matrix damage mode orthogonal to each other. The matrix layers between the laminae are resin-rich and can be modeled as isotropic materials; and hence, the interlamiar damage mode of delamination can thus be defined. As demonstrated in Figure 4.2, the intralaminar damage modes are modeled using the continuum damage model in this section and delamination is modeled using the CZM in ABAQUS/Standard. 4.4.2.2 Constitutive model To characterize the above damage modes, the damage tensor introduced in [7] is adopted here and is given by Eq (4.14). d1 0 0 D=diag{d1,d2,d3}= 0 d2 0 (4.14) o 0 d3 85 where the damage variable d1 is associated with fiber failure along the longitudinal direction, the damage variable d2 is associated with matrix failure (cracking or crushing) in the transverse direction, and the damage variable d3 is associated with delamination (In the implementation of current model, it is set to 0 because delamination is treated by CZM). The three damage modes are illustrated in Figure 4.7. (a) Fiber failure (b) In-ply matrix failure (c) Delamination Figure 4.7 The three damage modes in fiber-reinforced composites The actual stress and the effective stress introduced earlier now are tensors in three- dimensional. Let [c] be the actual stress tensor and [6] be the effective stress tensor as follows: 011 012 013 [6] = 0 22 0 23 (4.15) 'sym 0'33 and 511 5'12 5131' [61: 6'22 6'23 (4.16) Lsym 6'33 Using the transformation law proposed in [7], [6] can be obtained in terms of [o] and D as [6]=-;-[(1-D)"[a]+[o](1—D)"l (4.17) or explicitly 86 0'11 0'12( 1 + 1 ) 0'13( 1 + 1 ) 1-d1 2 1-d1 l—d2 2 l—d1 1—d3 [6]-- 3i “1% 1 + 1 ) (4.18) l—d2 2 l—dl 1—d3 ...... 1035 _ — 3 _ Using contracted notations, the above equation can be expressed as 6=Mo (4.19) where ._ . . . . . . T 6-{011,022,033,012»013,023} , (4.20) T 0={011.0221033312313923} (4,21) and the damage operator introduced earlier, which relates the actual stresses and the effective stresses, is now a tensor in three-dimensional and it takes the diagonal form 1 1 1 1 1 1 1 M: diagl —( ).-( )-( )} l-ld, ’1-c12’1-ar3 ’21-dl 1— d2 21—d1+1-d3 21- d2 1—d3 (4.22) Recalling that the effective stress is derived earlier using equilibrium equivalence (see Eq. (4.3)) in one-dimensional, the damaged constitutive equations in three-dimensional are now derived by the energy equivalence principle [7], which assumes that the elastic energy of the damaged material is in the same form as that of the undamaged material except the actual stresses are replaced by the effective stresses i.e. ‘P(6,dl) '—‘ “10(6) (4.23) where ‘P(o,d, ) = i—oTSdo (4.24) is the elastic energy in terms of the actual stresseso and the compliance matrix of the damaged material S d , and 87 ‘Fv-t- n- “I 1110(6) = £67806 (4.25) is the elastic energy in terms of the effective stresses 6 and the compliance matrix of the undamaged material S0. Inserting Eqs (4.24) and (4.25) into Eq. (4.23) and applying Eq.(4.19), we have Sd = MTSOM (4.26) Eq. (4.26) establishes the relationship between the compliance matrix of the damaged material Sd and that of the undamaged material S0. The relationship between the stiffness matrix of the damaged material and that of the undamaged material is derived by taking the inverse of Eq.(4.26) as follows. Cd = (Sd )‘1 = M‘1(SO)“M‘T = M"1C0M‘T (4.27) 4.4.2.3 Energy dissipation In case of elastic damaged material under isothermal conditions, the complementary free energy density takes the same form as the elastic energy in Eq. (4.24) in terms of the actual stresses. To ensure the thermodynamics irreversibility of the damage process imposed by the second law of dynamics, the rate of change of the complementary free energy minus the externally supplied work to the solid at constant strains must be non- negative: )? — (Fa 2 0 (4.28) Applying Eq. (4.24), it yields 67(51—a)+a—Wd, 20 (4.29) ad, 88 To ensure non-negative dissipation of mechanical energy, the above inequality must be satisfied for arbitrary stress rates. The expression in the parenthesis must be equal to zero, which yields 934:0 do This recovers the constitutive equation The inequality Eq. (4.29) becomes 3‘}! . . —d =Yd 20 ad] I I I where ,=av_=1,ras_d, 6d, 2 8d, are the thermodynamic forces conjugate to their (4.30) (4.31) (4.32) (4.33) respective damage variables d I (l=1,2,3), and they are associated with the three different damage modes. In the following, the thermodynamic forces will be shown to be non-negative. Taking the derivative of Eq (4.26), it gives 38,, 6M __ = MS __ ad, 0 ad, For a physical material So>O (4.34) (4.35) The damage operator M in Eq (4.22) is diagonal and positive definite, i.e. M>O Since 89 (4.36) 8M 1 —=———2 0 0110, (4.37) adl 2(1-41)2{ } 3M 1 = 020101, (4.38) adz 2(1—d2)2{ } and 9M=——1—2{0 0 2 0 1 1} (4.39) 6d3 2(1—d3) wehave 99120 1:123 (4.40) 3.1, Using Eqs (4.34), (4.35), (4.36) and (4.40), the following condition must hold fl 2 0 (4.41) 8d, Eq (4.33) gives Y, _>_ 0 (4.42) Therefore, the condition of positive evolution of damage variables (d1 20) is a sufficient condition for the proposed material model to fulfill the thermodynamics irreversibility imposed by the second law of thermodynamics. 4.4.2.4 Loading functions (damage initiation criteria) To simulate the progressive damage of laminated composite plates subjected to transverse loading, the loading functions F 1 (1:1, 2+, 2-, where 1 is associated with fiber damage, 2+ is associated with matrix cracking and 2- is associated with matrix crushing) are defined using the intralaminar damage initiation criteria developed by [21] as follows. For fiber failure 90 a 2 6'2 +a2 For transverse matrix cracking (6'22 2 0) 6' 6' 6' F2+ = J(—22 )2 +(——12 )2 +(———23 )2 (4.44) YT 512 Sm23 For matrix crushing (6'22 < 0) A 2A A A 1 a Y 0' 0' 0' F2_= _( 22)2+ C 22 _ 22 +( 12)2 (4.45) 4 512 4S122YC YC 512 where: 6'” - effective stress in the fiber direction (3'22 - effective stress in the transverse direction 612 - effective shear stress in the plane of fiber and transverse directions 6'23 - effective shear stress in the plane transverse and through-thickness plane 6'13 - effective shear stress in the plane of through-thickness and fiber directions X T - tensile strength in the fiber direction Y T - tensile strength in the transverse direction Y . . . . C - compresswe strength in the transverse direction S f - shear strength involving fiber failure $12 - shear strength in the fiber and transverse plane S m23 - shear strength for matrix cracking in the through-thickness direction 91 As can be seen from Eqs. (4.43), (4.44), and (4.45), the loading functions depend on the effective stresses and material properties. A value less than 1 of a loading function means no damage. A value greater than or equal to 1 of a loading function indicates that its associated damage mode is initiated. 4.4.2.5 Damage activation functions and consistency conditions The three damage activation functions f, (1:1, 2+, 2-) are defined using the loading functions as follows: f1 = F: - rim“ 5 0 (1:1. 2+. 2-) (4.46) where rim” 6 [1.r,f ] (1:1, 2+. 2-) (4.47) are the elastic domain thresholds associated with the damage modes 1:], 2+, 2- . They take an initial value of 1 before the material is damaged, and they increase with damage growth until reaching r,f (I=1, 2+, 2-) when the material is fully damaged. r1f is the elastic domain threshold associated with full damage in the fiber direction, rzf+ is the elastic domain threshold associated with full damage in transverse matrix cracking, and r2f_ is the elastic domain threshold associated with full damage in transverse matrix crushing. To ensure thermodynamics irreversibility, the damage activation functions Eq. (4.46) have to be always non-positive. While the damage activation function f, is negative, the material response is elastic. When a criterion is activated, f, = 0, the gradient of the loading function F, should be evaluated. If the gradient is not positive, the state is one of ' 92 unloading or neutral loading. If the gradient F, is positive, there is damage evolution; and the consistency condition has to be satisfied: f1 = F1 — f1max = 0 (4.48) The elastic domain thresholds take an initial value of 1 before material damage, and they increase monotonically with material damage. When the material is failed, they are set to be r,f . r,f (1:1, 2+, 2-) are the elastic domain thresholds corresponding to material failures. That is I‘lmax =min{max{l, max {Fl}}ir1f } (4.49) {=01max max _ - f r2+ —min max 1, max {F2+,F2_} ,r2+ (4.50) (=0Jmax rzl'llax =min max 1, max {F2_} ,rzji} (4.51) t=OJmax ,- where t is the loading time and rmax refers to the current loading time. Eqs. (4.49), (4.50), and (4.51) specify the condition that the elastic domain thresholds increase monotonically with material damage calculated through the loading functions F, . Eq. (4.50) takes account of the damage effect on matrix in tension due to matrix crushing. Two important characteristics of the model proposed here are that the elastic domain threshold values are a function of the damage variables and that the loading functions depend on the effective stress tensor. 93 4.4.2.6 Damage evolution and damage variable computation Similar to the one-dimensional case, once the material is damaged, the damage evolution is controlled by the following relation: :lf(r] max_ __1) d,- (4.52) ax(rlf - 1) In the above equation, r,f is calculated as follows: F of = ——2ECIS’ (4.53) (0'! ) wc where we is the crack band width, and the material properties associated with each damage mode are given in Table 4.1. Table 4.1 Material properties associated with each damage mode' in three- dimensional model Material Properties Damage Strength Young’s Moduli Fracture threshold toughness r,max 0'? El G ,F Fiber failure (I=l) rlmax X T E11 GIF Matrix cracking (I=2+) r ZinFaX YT 522 l G 5: Matrix crushing (I=2-) r2_ax YC 522 t G 2F_ To avoid the so-called snap-back effect in the stress-strain relation of the constitutive model, the following must hold: r,f 21 (4.54) which guarantees the strain energy released before damage initiation must be not more than the fracture toughness. 94 Since the finite element size must be not larger than the crack band width wc, the critical finite element size le can be obtained by combining Eqs. (4.53) and (4.54) as follows: F < 213,0, [6 S we _ (0f)2 (4.55) This critical finite element size (characteristic length) 18 is given as an input variable in ABAQUS/Standard user subroutine UMAT [1]. As shown later in the case study, this model does not eliminate mesh dependency completely; and it is recommended to align the mesh with the direction of crack propagation. 4.4.2.7 Viscous regularization of the damage variables Material softening and stiffness degradation in the elastic-softening model lead to severe convergence difficulties in implicit analysis programs, such as ABAQUS/Standard. In order to improve the convergence of the numerical algorithm, an artificial viscosity model originally proposed in [26] and used successfully for progressive damage analysis of composites in [12, 13] is implemented here as 4:” fun—db (4.56) where d)’ (I = 1,2+,2 —) are the “regularized” damage variables used in the calculation of the damaged elastic matrix, (I , (I =1,2+,2-) are the true damage variables, and ,u is the viscosity parameter representing the relaxation time to control the rate at which the regularized damage variables d }’ approach the true damage variables d, . Using viscous regularization with a small value of the viscosity parameter compared to the characteristic 95 time increment i.e./1 << At, enables the tangent stiffness matrix to be positive for small time increments during the softening regime, and thus, improves the rate of convergence of the model in the softening regime. As é—t —> oo, d 7 —> dv so that the solution of the ,u viscous system relaxes to that of the inviscid case. To update the “regularized” damage variables at time to +At , the above equation is discretized using the backward-Euler scheme in time as follows: At [1 dv : d +——dv 4.57 It0+At fl+Atlt0+At fl+AtltO ( ) Thus, ad; = A’ (4.58) 8d, fl+At Eq. (4.58) will be used in the calculation of the tangent constitutive tensor in the next section. 4.4.2.8 Solution algorithm and material tangent constitutive tensor The Newton-Raphson method is used to solve the system of nonlinear finite element equations. To ensure a fast convergence rate of the solution algorithm, the material tangent constitutive tensor must be calculated correctly, and it is computed as follows: V max A a.:_ d+z[a&8][ad, ad, ar, is] (4.59) I 3; - ad 1" ad, arlmax 86 as Using Eq (4.58) do At 3C d ad I dqmax as d ,u + At 7[ ad; 8134"“ aa 0 ( ) 96 where aCd is obtained by differentiating Eq (4.27), 8d, v max is computed from Eq. ad] 8r, max (4.52), and 2%— can be derived using the loading functions Eqs. (4.43), (4.44), and o (4.45). 4.4.2.9 Implementation of the material model The three-dimensional intralaminar continuum damage material model developed was implemented in ABAQUS non-linear finite element code using a user-subroutine UMAT and the flowchart is given in Figure 4.8. 4.4.2.10 Modeling of delamination Delamination is modeled using the three-dimensional cohesive zone model in ABAQUS with mixed-mode failure criteria, and the damage initiation and propagation criteria are given as follows. Damage initiation criterion [24] “my {if +[13—J2 =1 (4.61) K Tic 726 T36 where r, is the inter-laminar normal stress and 2'2 and T3 are the inter-laminar shear stresses, while If (i=1, 2, 3) are the inter-laminar strengths, and the symbol ( ) represents the Macaulay bracket defined as l- 0 Tl<0 . 97 L Enter UMAT ] , , Update strains and effective stresses Get material properties and "+1 n n+1 - - e = e + d8 compute stiffness matrix C0 ' 6714-1 =CO 811+] i Update thresholds Compute loading functions (rlmax )I‘H'l (Eqs 4.49_451) F] (Eqs. 4.43-4.45) d I = O No change to Cd Update damage variables d, (Eq. 4.52) Fully damaged: d, =1 Update C d (Eq. 4.27) Damage growing? No change to (rlmax )n-l-l > (rlmax )n d1 and Cd . . do Regularize damage variables a— : Cd (Eq. 4-60) (Eq. 4.57) 8 V f Compute damaged stiffness and tangent stiffness matrices Cd , g—a (Eqs. 4.27, 4.60) 8 Update actual stresses an+1 = Cd8n+l t [ Exit UMAT ] Figure 4.8 UMAT flowchart 98 The quadratic failure criterion Eq. (4.61) assumes that the compressive normal stresses do not affect the initiation of delamination. Delamination propagation criterion [25] {—01 ]+ {—011 J+[——G’” ]=1 (4.63) GIC GIIC Gmc where G,C(i = 1,11,!!!) are the critical energy release rates, which can be determined from experiments. G,(i = 1,11, 111) refer to the energy release rates. 4.5 Calibration of computational model To validate the model, a 90mm x 90mm square regular[0/ 9015 glass fiber reinforced polymer plate subjected to quasi-static displacement loading in [23] is simulated. The ply thickness is 1.07mm. Because of symmetry, only one quarter of the plate with a distributed load at a 6mm/sec rate in a central area of 0.5 mm x 0.5 mm is modeled (see Figure 4.9 (a)). Two zero-thickness cohesive layers of 25 mmx 25 mm are inserted at the center of the top and the bottom interfaces between the 0 and 90 plies (see Figure 4.9 (b)) to detect the delamination there. A mesh size of 0.5mm x 0.5 mm is used for the cohesive layers and the center of the plies, and a mesh size of 2mm x 2mm is used anywhere else. One element is used for each ply thickness. Due to identical mesh through thickness, shared node constraints are used to connect the cohesive layers and the plies. Incompatible mode continuum element C3D81 is used for the plies, and cohesive element COH3D8 is used for the cohesive layers. Clamped boundary conditions are imposed along the edges. The material properties are given in Table 4.2, where the fracture toughness associated with in-ply matrix failure is assumed to be the same as the Model I fracture toughness of the cohesive interface material properties from [17], and the 99 fracture toughness associated with in-ply fiber failure is assumed because of the lack of the data in literature. 1' " . -. p K‘ .“ .7 I >L .‘ ... .‘h ’v ~ ' :- y. "' .64“ ~ . ~ , 'j.’1 I .3“. ' ~. ‘ ‘ -. "11‘ “m » . O GSIVe a as ‘fr . . . 'Ja' .‘ \ 4 O a r l. 1 s 4» '0 b w‘ 4 x 1 - (a)Loading and boundary conditions (b) Delamination modeling Figure 4.9 Modeling of the quasi-static indentation test in [23] Table 4.2 Material properties used for three-dimensional progressive damage analysis GFRP ply material properties Ell I 522 I E33 G12 I 013 023 V12 V13 V23 (GPa) (GPa) (GPa) (GPa) (GPa) (GPa) 37.90 I 9.07 | 9.07 | 3.72 | 3.72 | 3.72 0.296 | 0.296 1 0.4 Strength Fracture toughness (GPa) (J/mz) - , 1 X7 IYTI YC I 512 I Sf I523I 61F I 02’: I Of. 900 | 74 | 237 | 64 | 190 | 64 l 2880 | 240 | 240 Cohesive layer material properties Stiffness (Calculated) Strength Fracture toughness (MP3) (MP3) ( J lmz) k‘ I k2 I k3 If I 15 I I; GICI GIIC I 01116 24x106 I 34,5,(106 I 34.5x106 24.0 I 42.7 I 42.7 240 I 640 I 640 The numerical results correlate well with the experimental results. The load- displacement curve comparison is given in Figure 4.10. The comparison of delamination shape and size at the bottom cohesive layer is given in Figure 4.11. Interestingly, if only the delamination failure criteria for the cohesive layers are activated and the intralaminar 100 damage is not allowed, the delamination will grow into an ellipse as illustrated in Figure 4.12. From the comparison of the results with intralaminar damage criteria activated and deactivated, we can conclude that the intralaminar damage, matrix cracking, along the fiber direction in the bottom ply promotes the delamination to grow to the peanut shape obtained from experiment. 3500. . . - . . a - r . . . . 3000. . _ ,3"... 1 2500- " . ... 2000. _ . 3. “- 1500.» . -' Experimental ---- NumerlcaltCurrent) 1000.» - 500- 1 O, i i . 1 i i i i i i I i i . 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Displacementlmm) Figure 4.10 Load-displacement curve comparison 101 (a) Simulated delamination (39.2mm) (b) Matrix cracking at bottom ply (c) Experimental delamination Figure 4.11 Delamination comparison between numerical and experimental results in [23] Figure 4.12 Delamination with intralaminar damage deactivated 102 4.6 Case study-impact induced damage modes This is a qualitative study of the delamination patterns and fiber—bridging damage mode in two-ply skew-angled [0/01] (01:90, 75, 60, 45, 30, 15, and 10 degree) laminated plates (see Figure 4.13) subjected to transverse low velocity impact. Delamination patterns in laminated plates subjected to transverse low velocity impact are shown to be independent of material properties in the original experimental work by [2]. Material properties used here are the same as for the model calibration given in Table 4.2. Top 0 ply Cohesive layer Bottom (1 ply Figure 4.13 Modeling of [Old] laminate 4.6.1 Model I In Model I, the whole 90mm x 90 mm square plate (see Figure 4.14 (a) ) is modeled to study the skew angled laminates. Clamped boundary conditions are used, and uniform transverse displacement at a 6mm/sec loading rate is applied to the 9 nodes at the center. A 50mm x 50mm zero-thickness cohesive layer is inserted into the center of the interface between the t0p and bottom plies (see Figure 4.14 (b)). The mesh used is shown in Figure 4.14 (c), a mesh size of 0.5 mm x 0.5 mm is used for the cohesive layer and the center of the plies, and a mesh size of 2mm x 2mm is used everywhere else for the plies. The ply 103 thickness is 0.5 mm and one element is used for each ply thickness. Again the elements used are C3D81 for the plies, and COH3D8 for the cohesive layer. The total element number in the model is 38,800. (a) Loading and boundary conditions 0 ply ": my Cohesive layer 6 (b) Cohesive layer (c) Mesh Figure 4.14 Model I with the same mesh through thickness 104 4.6.1.1 Results and discussions In Figure 4.15, the delamination patterns obtained using Model I for [0/01] (01:90, 75, 60, 45, 30, and 15 degree) laminated plates are compared with experimental results from [2]. For large angle (a 2: 30 degree) laminates, the correlation between experimental and numerical results is very good. However, the numerical delamination shape for small angle [0/ 15] laminate is not observed in experiment. This numerical artifact is due to mesh sensitivity. The damage modes in the small angle laminates will be investigated using Model 11 introduced next. In large angle (a 2 30 degree) laminates, delamination and matrix cracking in the bottom plies are the dominant damage modes, and the delamination in the interface layers and the matrix cracking in the bottom layers are shown in Figure 4.16. Similar to the result in the cross-ply laminate, the matrix cracking in the bottom plies promotes the delamination to grow into the peanut shapes along the fiber directions in the bottom plies. 105 [0/90] a [0/60] ' Figure 4.15 Delamination pattern comparison between experimental (center) from [2] and numerical (left and right) using Model I [O/90] ' [0/75] ‘ [0/60] [0/30] I . l [0/901 l [0/75] I [0/60] I [0/45] . (a) Delamination Matrix cracking at bottom layers .- [0/30] [0/15] I Numerical artifact i because of mesh sensitivity Figure 4.16 Damage modes in large-angle [0/ a] laminates (a>=30) 106 4.6.2 Model 11 To investigate the damage modes in small angle laminates, Model 11 as illustrated in Figure 4.17 (a) with a circular [0/ 15] plate of 90mm in diameter is employed. Clamped boundary conditions are used, and uniform transverse displacement at a 6mm/sec loading rate is applied to the 9 nodes at the center. A zero-thickness cohesive layer of 50mm in diameter (see Figure 4.17 (b)) is inserted into the center of the interface between the top and bottom plies. The ply thickness used is 0.5mm, and one element is used to model each ply in thickness direction. To model the matrix cracking along the fiber direction in each ply, a mesh size of 0.5mm x 0.5mm is used for the rectangular area of 32mm x 50mm in the center of the plies and the cohesive layer. The mesh direction is aligned with the fiber direction in each ply as shown in Figure 4.17 (c) and (d). Automatic mesh scheme is used everywhere else. The elements used are C3D81 for the plies and COH3D8 for the cohesive layers. The total element number in the model is 31,553. Surface-based tie constraints are used to tie the cohesive layer with the plies. 4.6.2.1 Results and discussions As explained in the development of the damage model, the crack band model does not eliminate mesh dependency completely; and it is recommended to align the mesh with the direction of crack propagation. Using Model II with mesh aligned with ply fiber orientation enables the correct prediction of damage mode in small angle laminates as illustrated in Figure 4.18 because cracks propagate along the fiber direction in fiber- reinforced composites. As shown in Figure 4.18, the delamination using Model II is negligible. This shows that the crack band model can overcome mesh sensitivity due to size, but it can not avoid the mesh sensitivity due to direction. The top views of damage 107 evolution in a [0/15] laminate given in Figure 4.19 show that matrix cracking first initiates along the fiber direction in the center line of the bottom ply due to tension, and then matrix cracking along the matrix cracking line in the top ply follows. This happens because the two plies are bonded together due to negligible delamination. This damage mode is called fiber-bridging, and it can be used to combat delamination using laminates with small angle between adjacent plies. The fiber-bridging damage mode is illustrated using a [0/10] laminate in Figure 4.20. Detailed discussion of its advantages can be found in the experimental work by Pilchak et al. [3]. I. Cohesive layer -, . 0 ply ‘15 ply (a) Loading and boundary conditions (b) Cohesive layer (c) Top layer ((1) Bottom layer Figure 4.17 Model 11 with mesh aligned with ply fiber orientation 108 [0/15] 7 Delamination Delamination Matrix cracking Matrix cracking at bottom 15 ply at bottom 15 ply (a) Model I (b) Experiment (0) Model II Figure 4.18 Result comparison of [0/15] laminate using Model I (left) and Model H (right) with experiment (center); Pictures in the top show the delamination and pictures in the bottom show the matrix cracking in the bottom ply. (b) Delamination (c) Bottom ply (15) with Matrix cracking Figure 4.19 Damage evolution in [0/15] (Top view) 109 (a) Matrix cracking in top ply(0) (b) Delamination (0) Matrix cracking in bottom ply(10) Fibers in t0p layer 0 - N (0 ply) form bridges Matrix cracking in 0 \ bottom layer (10 ply) creates a river (d) Fiber-bridging diagram e) Fiber-bridging from experiment [3] Figure 4.20 Fiber bridging in a [0/10] laminate 4.7 Stacking sequence effect on the delamination area Figure 4.21 shows the comparison of predicted normalized delamination areas between numerical simulation and bending stiffness mismatch theory [2] in two-ply [0/01] laminates with respect to the change of stacking angle 01. In Figure 4.21, all the 110 delamination areas are normalized with respect to the delamaintion area of [0/90] laminate. For 01 >15, the predicted delamination areas from both the numerical simulation and the bending stiffness mismatch theory increase with the increase of the angle 01. However, no delamination is predicted in the numerical simulation for 01 <15 because of the fiber-bridging effect, which is in line with the experimental results in [3]. 1 a" +Numerical Results I __ __ .— .— -— / — —- Bending Stiffness / § 0-8 7” Mismatch Theroy / 8 < g 6 0.6 —- w ._ (u 25:5 : E52 0.4_A A A - O E / Z .1..“ / (D / D 0.2 T‘F ’ I V - 4— 4—4 72L __ O ’ ’ l i l l 1 0 10 15 30 45 60 75 90 a(Degree) Figure 4.21 The comparison of predicted normalized delamination areas between numerical simulation and mismatch theory [2] 4.8 Conclusions of three-dimensional damage analysis of laminates o A three-dimensional continuum damage model for the prediction of the initiation and evolution of intra-laminar damage mechanisms is developed by integrating the stress-based failure criteria for damage initiation, fracture-mechanics based criteria for damage propagation, and the crack band model for alleviating mesh sensitivity. Thermodynamics irreversibility of this model is ensured by the positive evolution of damage variables. 111 The developed intralaminar damage model combined with the cohesive zone model in ABAQUS enables the correct prediction of damage modes (i.e., delamination and fiber-bridging) in laminates subjected to transverse static loadings. Damage in laminates with large angles (0t>=30) between adjacent plies is dominated by delamination, and the matrix cracking along fibers in the bottom ply promotes the delamination to grow into peanut shapes. Damage in laminates with small angles (a<20) between adjacent plies is dominated by fiber bridging, which is a preferred damage mode over delamination. To predict the fiber bridging damage mode in small angle laminates, mesh aligned along fibers in each ply should be used. Severe convergence difficulties are observed in the three-dimensional progressive damage analysis using ABAQUS/Standard. The analysis is often terminated before completion of the analysis job. The results presented are recorded at different loading times in the loading history. A 2-ply laminate modeled with the cohesive zone model only inserted in the center results a total element number of between 30,000 and 40,000 in the model. Typically, a job needs about one week to complete using a computer on a Linux machine with a Dual Intel Xeon E5345 2.33GHz 64bit processor and 8GB RAM. Therefore, this model must be used selectively for computation efficiency. 112 REFERENCES . ABAQUS/Standard User’s Manual, online ver. 6.8 (2008), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. D. Liu, Impact-Induced Delamination—A View of Bending Stiffness Mismatching, Journal of Composite Materials, Vol. 22, No. 7(1988), 674-692 AL. Pilchak, T. Uchiyama and D. Liu, Low-velocity impact response of small-angle laminated composites, AIAA Journal, vol. 44, No.12 (2006) 3080-3087. J. Lemaitre, A course on damage mechanics, Berlin ; New York : Springer, 1996, 2nd rev. L.M. Kachanov, On the creep fracture time, Izvestiya Akademii, Nauk USSR Otd. Tech. 8 (1958), pp. 26—31 (in Russian). Y. N. Rabotnov, Creep rupture. Proc. XII, Int. Cong. on Applied Mechanics, StanfordSpringer, 1968. J .P. Cordebois and F. Sidoroff, Damage induced elastic anisotropy. In: J .P. Boehler and Martinus Nijhoff, Editors, Mechanical Behavior of Anisotropic Solids, The Hague (1982) 761—774. ' A. Matzenmiller, J. Lubliner and RL. Taylor, A constitutive model for anisotropic damage in fiber-composites, Mech Mater 20 (1995) 125 - 52. M. Zako, Y. Uetsuji, and T. Kurashiki, Finite element analysis of damaged woven fabric composite materials, Composites Science and Technology 63 (2003) 507—516 10. P. Maimi, PP Camanho, J Mayugo and CG Davila, A thermodynamically consistent 11. 12. 13. 14. damage model for advanced composites, Tech. Rep. NASA/1‘ M-2006-214282, 2006. P. Maimi, PP Camanho, J Mayugo and CG Davila, A continuum damage model for composite laminates — part I: constitutive model. Mech Mater 39 (2007) 897—908. P. Maimi, PP Camanho, J Mayugo and CG Davila, A continuum damage model for composite laminates — part II: computational implementation and validation.Mech Mater 39 (2007) 909—19. 1. Lapczyk and J. A. Hurtado, Progressive damage modeling in fiber-reinforced materials, Composites: Part A 38 (2007) 2333—2341 F. Aymerich, F. Dore and P. Priolo, Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements, Composites Science and Technology 68 (2008) 2383—2390 113 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. M. Nishikawa, T. Okabe and N. Takeda, Numerical simulation of interlaminar damage propagation in CFRP cross—ply laminates under transverse loading, International Journal of Solids and Structures 44 (2007) 3101—31 13 Z. Zou, S.R. Reid, S. Li, and PD. Soden, Application of a delamination model to laminated composite structures, Composite Structures 56 (2002) 375—3 89 N. Hu, Y. Zemba, T. Okabe, C. Yan, H. Fukunaga, and AM. Elmarakbi, A new cohesive model for simulating delamination propagation in composite laminates under transverse loads, Mechanics of Materials 40 (2008) 920—935 C.F. Li, N. Hu, Y.J. Yina, H. Sekinec, and H. Fukunagac, Low-velocity impact- induced damage of continuous fiber-reinforced composite laminates. Part I. An FEM numerical model, Composites: Part A 33 (2002) 1055-1062 C.S. Lopes, P.P. Camanho, Z. Gtirdal, P. Maimi’, and E.V. Gonzalez, Low-velocity impact damage on dispersed stacking sequence laminates. Part 11: Numerical simulations, Composites Science and Technology 69 (2009) 937-947 Z.P. Bazant and J. Planas, Fracture and size effect in concrete and other quasi-brittle materials, CRC Press, 1998 Z.P. Bazant, and BH. Oh, Crack band theory for fracture of concrete, Materials and structures 1983 (16) 155—177 J. P. Hou, N. Petrinic, C. Ruiz, and S. R. Hallett, Prediction of impact damage in composite plates, Composites Science and Technology 60 (2000) 273—281 S. Kamiya, H. Sekine and Y. Yagishita, Computational simulation of interlaminar crack extension in angle-ply laminates due to transverse loading. Journal of Composite Materials 32 (1998), 744—765. W.,Cui, M. R. Wisnom, and M. Jones ”A Comparison of Failure Criteria to Predict Delamination of Unidirectional Glass/Epoxy Specimens Waisted Through the Thickness.” Composites 23(3), .1992,158-66. E.M. Wu and Jr C. Reuter, Crack extension in fiberglass reinforced plastics, Report no. 275, T&AM, University of Illinois, 1965 114 Chapter 5 Three-dimensional progressive damage analysis of woven composites subjected to impact loading using shell-solid coupling technique 5.1 Abstract The objective of this chapter is to evaluate the dynamic response of laminated and woven composite plates subjected to transverse low-velocity impact. The shell-solid coupling technique is implemented to model the complex geometry of woven composites because of the numerical accuracy of the three-dimensional solid elements and the computational efficiency of the shell elements. The studies use the homogenization method developed in Chapter 2 to compute the effective material properties for the global shell elements. To model delamination and intra-laminar damage modes in the local region, the cohesive zone model (CZM) presented in Chapter 3 and the continuum damage mechanics (CDM) model developed in Chapter 4 are employed. 5.2 Introduction Compared with the progressive damage analysis for laminated composites, the progressive damage analysis for woven composites has not been well investigated. Currently, studies of delamination of woven composites using the cohesive zone model (CZM) and continuum damage mechanics (CDM) model have been limited to the sublarninate model with each layer modeled as homogeneous and orthotropic material. Using the sublarninate model, lannucci et. al. [1-3] studied the damage modes in plain woven composites subjected to impact loading. In their study, the laminated plain woven composites were modeled with two homogeneous layers and a resin-rich interface in between. The damage modes of homogeneous layers were modeled with plane stress shell elements, and delamination on the interface was modeled with cohesive elements. 115 The sublarninate model using shell elements is efficient. However, the information between plies of each woven layer can not be obtained because homogeneous material properties are assumed for each layer. In order to obtain the damage information within each woven layer, full three-dimensional modeling of woven composites may be necessary. However, full three-dimensional progressive damage analysis of woven composites poses the following challenges: 1. The finite element modeling of woven composites is much more complex than the modeling of laminates due to the irregular matrix pockets and the undulation regions present in woven composites. Different elements and meshing schemes have to be used for the laminated region, undulation region, and pure matrix pocket region. While the cohesive zone model (CZM) may be useful for modeling the fracture interfaces between adjacent layers of laminates, convergence problems due to material softening make implicit solvers, such as ABAQUS/Standard, impractical for progressive analysis of woven composites. Furthermore, inserting zero- thickness cohesive layers into different locations of woven composites is also a challenge. The Continuum Damage Mechanics (CDM) model may be a viable solution for woven composites. However, it has the following potential problems: 0 Similar to the CZM, the material softening behavior in CDM model is known for possible lead to severe convergence difficulties in implicit solvers such as ABAQUS/Standard. In implicit solvers, iterative algorithms such as the Newton-Raphson method are used to solve the non-linear system of 116 equilibrium equations. With material softening and stiffiiess degradation in a model, there is a tendency for damage to accumulate in small elements. The degradation of stiffness due to accumulated damage may cause the tangent stiffness matrix to be non-positive definite, leading to severe convergence difficulties. 0 In a CDM model with fracture mechanics based damage propagation criterion, if the finite element size is larger than the smearing crack band width, the dissipated energy (strain energy release rate) may become greater than the fi‘acture toughness of the material when stresses reach the strength of the material. The stress will then drop to zero immediately. The finite element analysis will become unstable leading to the so-called snapback issue. To avoid the snapback issue in models with fracture mechanics based criteria for crack propagation, the finite element size must be small enough, typically smaller than 0.5mm for fiber-reinforced polymer matrix composites. This imposes high demands on computation resources and efforts. 4. The experimental methods for obtaining some of the material properties including fracture toughness and strength employed in the damage initiation and propagation criteria in composites are still under development. This poses another challenge for the progressive damage analysis of composites using the fi'acture mechanics based criteria as those experienced in. CZM and CDM models. Published data of these material properties are very limited. Some of the material properties used in this study are from the literature and some others are assumed 117 as explained in the three-dimensional static progressive damage analysis in Chapter 4. 5.3 Proposed Methodology To address the above challenges of progressive damage analysis of woven composites, the organization of dynamic progressive damage analysis shown in Figure 5.1 is adopted in this study. To model the complex geometry of woven composites, a shell-solid coupling technique is used. The computed effective material properties using the homogenization method developed in Chapter 2 are assigned to the shell sections of the global region. Regarding the modeling of the damage modes in the local region, delamination is handled with the cohesive zone model presented in Chapter 3, and the intralaminar damage is tackled with the continuum damage mechanics model developed in Chapter 4. The detailed discussions of the modeling techniques are followed. 5.3.1 Shell-solid coupling technique In light of the complexity of the modeling of woven composites, the shell-solid coupling technique available in ABAQUS will be employed. The shell-solid coupling technique combines the numerical accuracy of three-dimensional solid elements and the computational efficiency of shell elements. Krueger et al. [4-7] have applied shell-solid coupling technique successfully in the investigation of delamination in laminated composite structures. Good correlation with experiment and simulation based on full three-dimensional model has been achieved. In the current study of composites subjected to impact loading using the shell-solid coupling technique, each ply in the local areas of interest in the composite specimens is modeled with continuum elements, interfaces 118 between plies are modeled with cohesive elements; and the rest of the composite specimens are modeled with shell elements. Modeling Dynamic Progressive Damage Analysis in ABAQUS/Explicit | f i Geometric modeling using I Material modeling I shell-solid coupling I r4“: """"" l l Model global Modeling of Homogenized Continuum Damage elastic regions Local regions material property Mechanics Model with shell with damage for global shell for local solid elements I elements (Chap. 2) elements (Chap.4) * . £ . a Cohesive Zone Model I Solid elements I Cohesrve elements I‘ _ _ _ for local cohesive for plies for interfaces elements (Chap. 3) Figure 5.1 Organization of dynamic progressive damage analysis 5.3.2 Material modeling The three-dimensional continuum damage model developed in Chapter 4 for the prediction of the initiation and evolution of intra-laminar damage mechanisms will be implemented in ABAQUS/Explicit [8] using a user subroutine VUMAT. The VUMAT will be used to control the initiation and evolution of damages within each ply. The regions where no damage is expected will be modeled with shell elements, and the material properties computed using the homogenization technique developed in Chapter 2 will be assigned to these shell elements to achieve high computational efficiency. 119 5.3.3 Solution technique To overcome the convergence problems in implicit solvers caused by material softening present in the CZM and CDM, ABAQUS/Explicit solver will be used for the dynamic analysis in this study. 5.3.4 Cohesive layer modeling Damaged interfaces between different layers in laminated and plain woven composites are modeled with zero-thickness cohesive layers. Damaged interfaces with undulations in plain woven and Q3D woven composites are difficult to model with zero- thickness cohesive layers; and thus, they are modeled with a thickness of 0.1% of the ply thickness through PYTHON scripts in ABAQUS/CAB. In the delamination model using cohesive elements with constitutive response defined in terms of traction-separation laws as discussed in Chapter 3, the cohesive layer represents an infinitesimally thin layer of adhesive, typically modeled with geometric thickness close or equal to zero. The geometric thickness does not affect the solutions because the constitutive thickness of 1.0 is used for traction-separation response in ABAQUS to ensure that the nominal strain is equal to the separation (relative displacement). However, typically zero thickness or less than 3.0% of ply thickness is used by most researchers, e. g. [9] in modeling the cohesive layers. 5.3.5 Implementation of the material model The three-dimensional intra-laminar continuum damage model developed in Chapter 4 was implemented in ABAQUS/Explicit non-linear finite element code using a user- subroutine VUMAT and the flowchart is given in Figure 5.2. The notations and the 120 equations used here are the same as those shown in the UMAT flowchart given in Figure 4.8. 5.3.6 Modeling of delamination The interlaminar damage (delamination) is modeled using the three-dimensional cohesive zone model in ABAQUS with mixed-mode failure criteria. The damage initiation and propagation criteria are the same as in Chapter 4. 5.4 Model calibration 5.4.1 Material model calibration A uniaxial tension test is carried out using a 0.25m x 0.25m x 0.25mm cube modeled with one continuum three-dimensional C3D8R element with reduced integration to verify the material model implemented in VUMAT. As illustrated in Figure 5.3, uniform static displacement loading along the vertical direction is prescribed at the nodes of the top surface of the cube to simulate homogeneous deformation and stress conditions. Zero displacements along the vertical direction are assigned to the nodes on the bottom surface. To prevent rigid body motion modes, zero displacements in the other two directions are also imposed on one of the corner nodes at the bottom surface. This model is analyzed using the VUMAT in ABAQUS/Explicit and in ABAQUS/Standard (an implicit solver) using the UMAT developed earlier in Chapter 4. The resulting force-displacement curves are compared in Figure 5 .4. To make sure the energy in VUMAT is updated correctly, the internal energy is plotted and compared in Figure 5.5. As can be seen, the results fi'om the UMAT in ABAQUS/Standard (implicit solver) correlate very well with those of the VUMAT in ABAQUS/Explicit. 121 1? [ Enter VUMAT J , , Update strains and effective stresses Get material properties and "+1 n n+1 - - a = a + d3 compute stiffness matrix C0 6n+l =CO 8n+l 1, Update thresholds Compute loading functions (rlmax)n+1 (Eqs. 449451) F, (Eqs. 4.43-4.45) d, = 0 No change to Cd Update damage variables d, (Eq. 4.52) Fully damaged: d, =1 Update Cd (Eq. 4.27) Damage growing? No change to (rlmax)n+l >(r1rnax)n d, and Cd Y , Update actual stresses Update damaged matrix "+1 "+1 Cd (Eqs. 4.27) - 5 = (348 v I Exit VUMAT I Figure 5.2 VUMAT flowchart 122 Figure 5.3 Static tension test on a single element 1.2 0.84 5 — Explicit § Implicit o u. 0.45 0.0 , , , , , - , - , 0.0 1.0 2.0 3.0 4.0 5.0 6.0 [x1 .E-S] Displacement (m) Figure 5.4 Force-displacement comparison 123 [x1.E-6] 7.0 ~ 6.0r 5.0 T — Explicit 4.0 - ..... Implicit 3.0- Energy (J) - 2.0- 1.0- 0.0. ----- fii 0.0 110 ' 2I0 ' 310 ' 4J0 ’ 510 ' 6.0[x1.E-6] Displacement (m) Figure 5.5 Internal energy comparison 5.4.2 Dynamic modeling calibration To ensure that the impact model including contact algorithm and dynamic solver are set up correctly, a transverse impact test is performed on an aluminum plate and compared with the experiment [10]. As illustrated in Figure 5.6, a 76.2mmx 76.2mm x 2.06mm (length x width x thickness) plate made of A1 6061 T6 is modeled as an elastic- plastic material using the J ohnson-Cook plasticity model available in ABAQUS/Explicit. In this model, the whole plate is modeled to compare the displacement contour later on. The material properties used are as follows: Young’s Modulus E=68.9GPa, density p =2720 kg/ m3 , Poisson’s ratio v =0.33, and the J ohnson-Cook Coefficients from [11] are used for the J ohnson-Cook plasticity model. The impactor with a hemisphere of 12.7mm in diameter is modeled as a rigid analytical surface associated with a point-wise mass of 5 .2 kg, which are the same as in the experiments. An initial velocity of 1.62m/s in the impact direction is prescribed to the 124 impactor, simulating the impact velocity measured during the tests. Additionally, a gravitational force is applied to the impactor to take account of gravity. The general contact algorithm in ABAQUS/Explicit, which uses a penalty enforcement contact method, is employed to simulate the contact between the impactor and the specimen. The friction coefficient of 0.1 is assumed and introduced between the impactor surface and the specimen. The plate is modeled with elements of type C3D8R, and clamped boundary conditions are used. A fine mesh size of 1.27mm x 1.27mm is used for the area under impact, and coarse mesh size of 6.35m x 6.35mm is used for the rest of the plate. The model is analyzed in ABAQUS/Explicit, and comparisons with experiments of the load history, load-displacement curve, and displacement contour are given in Figure 5.7, 5.8, and 5.9, respectively. Reasonable overall correlation between numerical simulation and experimental result has been achieved. The unloading curves match well between the numerical and the experimental results. The over-prediction of the peak force and displacement may be due to the material properties used. As explained earlier, material properties fi'om Ref. [11] were used in the simulation due to the lack of the tested material properties available. Figure 5.6 Transverse impact of aluminum plate model 125 4500 o i V i 1 fl 1 4000. - 3500. 3000. p I g 2500. a 2000. I 1500. ,- _ Experimental — Numerical 1000. 500. 0' . 1 .¥ I L l A l . l i 1 0.000 0.001 0.002 0.003 0.004 0.005 0.006 Time(second) Figure 5.7 Comparison of force history with experiment [10] 4050 ‘ I ' T ' I r I ‘ I " T ‘ l ‘— 4.00 _ Experimental —-- Numerical r—‘Wr 3.50 1" U! o , . Foroe(kN) l" o o V 1.00- 050- A, 0.00 . . . . - . . . . . . L 0.0 0.5 1.0 1.5 ‘ 21.0 ' 2.5 3.0 3.5 ‘ 4.0 Dlsplacement(mm) Figure 5.8 Comparison of force-displacement curve with experiment [10] 126 (a) Numerical (b) Experimental[10] Figure 5.9 Comparison of displacement contour with experiment [10] 5.4.3 Shell-solid coupling model calibration To ensure that the shell-solid coupling technique is applied properly, the quasi-static indentation test on the 90mm x 90mm square [0/90]s glass fiber reinforced polymer plate [12] studied in Section 4.5 is employed again here for the shell-solid coupling model calibration purpose. Due to material and geometric symmetry, only a quarter plate of 45mm x 45mm is modeled using the shell-solid coupling technique as shown in Figure 5.10. In the shell-solid model, the square of 25mm x 25mm in the center of the plate where damage is expected to occur is modeled with three-dimensional elements and the rest of the model is modeled with plain stress shell elements 84R. 127 Simply supported along the edges Shelll-solid / Global elastic region 001113 111? modeled with shell elements constraint S4R Plies with damage at local regions modeled With 5059 Interfaces with delamination elements C3D8R at local regions modeled S etries along the 331111 $538M: elements symmetric planes Figure 5.10 Shell-solid coupling model of [0/90]s plate for static analysis In the local three-dimensional model, solid elements of type C3D8R are used for the plies and the three-dimensional cohesive elements COH3D8 are used for the interfaces between plies. Again, two cohesive layers are inserted into the interfaces between the 0° and 90° plies of the local three-dimensional model. Enhanced hourglass control based on the assumed enhanced strain method for solid elements is used for the element type C3D8R. A mesh size of 0.5mm x 0.5mm is used for the solid section and mesh size of 2mm x 2mm is used for the shell section. Simply-supported boundary conditions are applied along the edges of the shell, and symmetries are imposed on the associated symmetric boundaries of the quarter plate. Displacement loading is imposed on the center 4 nodes of the model, similar to that in Section 4.5. The material properties computed using the homogenization technique developed in Chapter 2 are given in Table 5.1. They are assigned to the global shell section. The material properties to be used for the local solid and cohesive sections are given in Table 128 3 5.2. Table 5.2 is the same as Table 4.2.Besides, density p = 1850 kg/ m and p = 1210 kg/ m3 are used for the composite plies and cohesive layers, respectively. Table 5.1 Homogenized material properties for global shell section of [0/90]s V13 = V23 I 1511 = E22 E33 012 613 = G23 V12 (GPa) (GPa) (GPa) (GPa) | 23.58 | 9.08 | 4.14 | 3.75 | 0.09 | 0.36 | The model is analyzed using ABAQUS/Standard. The results fiom the shell-solid model correlate well with those fiom experiment [12] and full three-dimensional model (see section 4.5). The comparison of load-displacement curves is given in Figure 5.11. As can be seen from Figure 5.11, the full three-dimensional model over-predicts the force when displacement is greater than 2.3mm. This may be due to the use of element C3D81 with full integration. The shell-solid model (using element C3D8R with reduced integration) denoted as Shell/3D in the plot under-predicts the force. The difference increases with the growth of the displacement. This may be because the shell-solid coupling constraint becomes troublesome as the damage (matrix cracking and delamination) gets close to the coupling boundary. The comparison of the forces obtained at the delamination length of 40mm is given in Table 5.3. The comparison of delamination shape and size at the bottom cohesive layer is given in Figure 5.12. No delamination at the top cohesive layer was reported in the experiment [12]. The damage evolution and the predicted unphysical delamination at the top cohesive layer are discussed next. 129 Table 5.2 Material properties used for local solid and cohesive sections GFRP ply material properties E11 522 E33 012 013 023 V12 V13 V23 (GPa) (GPa) (GPa) (GPa) (GPa) (GPa) 37.90] 9.07 | 9.07 | 3.72 | 3.72 | 3.72 0.296 | 0.296 | 0.4 Strength Fracture toughness (GPa) (flmz) XT I YT I YC I S12 I Sf |S23 GIF I 02]: I 65: 900 | 74 | 237 | 64 | 190 | 64 2880 | 240 | 240 Cohesive layer material properties Stiffness (Calculated) Strength Fracture toughness (MPa) , (MPa) (J/mz) 1‘1 I k2 I 1‘3 1f I a; I a; GICI GHC IGIIICI 24x106l34.5x106I34.5x106 24.0] 427 I427 240| 640 I640 3500. . 4 - 3000.. ----- Full 30 — Experiment ......... 2500,. — Shell/3D .......... , 22000. ....... , O J 2 11°. 1500. 1000. 500. o. .1 . - . . ._ . i 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Displacement (mm) Figure 5.11 Force-displacement curve comparison Table 5.3 Comparison of forces obtained at delamination length of 40mm Experiment Numerical [12] Full three- Shell/three-dimensional dimensional Force (kN) 2.72 2.85 2.46 Deviation (%) N/A 4.8 -9.6 130 (3) Predicted delamination at bottom 90/0 interface (b) Predicted matrix cracking at bottom 0 ply . w...“ I W4 i . ‘7 a. . i40mm { d (c) Delamination at bottom 90/0 interface from experiment [12] Figure 5.12 Comparison of damage The damage evolution diagrams are given in Figure 5.13, while the interlaminar stresses of the top 0/90 interface and the bottom 90/0 interface are given in Figure 5.14 and Figure 5.15, respectively. With the comparisons of Figure 5.13 (a) and Figure 5.14 for the top 0/90 interface and Figure 5.13 (b) and Figure 5.15 for the bottom 90/0 interface, it can be seen that the maximum interlaminar shear stresses occur at the delamination fronts (fracture process zones). The matrix cracking at the centerline in the bottom 0° - ply and the peanut-shaped delamination at the bottom 90/0 interface correlate very well with the experimental results given in Figure 5.12 (0). However, the delamination at the top 0/90 interface as shown in Figure 5.13 (a) is much larger than that observed in experiments. This may be due to the inaccuracy of the failure criterion used for the delamination prediction in this study. 131 (a) Delamination at top 0/90 interface (b) Delamination at bottom 90/0 interface (c) Mam'e'r‘aiagg Ofbottom 05.31;--- _. . Figure 5.13 Damage evolution (only a quarter of specimen shown) 132 5, $33 (Avg: 75%) -5.968e+07 -1.3190+08 -2.041e+08 -2.764e+08 -7.820e+08 -8. 5428+08 (a) Interlaminar normal stress 5, 513 S, 523 (Avg: 75%) (Ava: 75%) +3.1809+07 +3.327e+07 +2.746e+07 +2.964e+07 +2.312e+07 +2-6028+07 +1.878e+07 +2-2409+07 +1.445e+07 +1.877e+07 +1.011e+07 +1.515e+07 +5.770e-l-06 +1.1526+07 +1.432e-I-06 +7.900e+06 - .906e +4.27Ge-i-06 -7.244e+06 +6.521e+05 -1.158e+07 -2.972e+06 -1.592e+07 -6.596e+06 -2.026e+07 -1.022e+07 (b) Interlaminar shear stresses Figure 5.14 Stresses at top 0/90 interface 133 S, 533 (Avg: 75%) +1.133e+07 -1.603e+07 -4.340e+07 -7.077e+07 -9.813e+07 -1.255e+08 - 1.529e+08 -1.8029+OB -2.076e+08 -2.350e+08 -2.623e+08 -2.897o+08 -3.171e+08 (a) Interlaminar normal stress S, 513 5' 523 (Avg: 75%) (Avg: 75%) +3.3480-l-07 +2.785e+07 - +3.01%:3; *i'éiézigi hiilammum + 5399 + . 4.2.3599“); +1.32%”, 3 ear stresses +2.030e+07 +8.348e+06 +1.700e+o7 +3.474e+06 +1.370e+07 -1.400e+06 +1.041o+o7 -6.Z75e+06 +7.110e+06 -1 115e+07 +3.813e+06 -1 602e+07 +5.163e+05 -2.090¢+07 -2.7800+06 -2.577e+07 -6.077e+06 -3.065e+07 (b) Interlaminar shear stresses Figure 5.15 Stresses at bottom 90/0 interface 5.5 Case study of impact- induced damage 5.5.1 Comparison of laminate, plain weave, and Q3D weave specimens In this section, the laminate L[0/9O]6, plain weave 2D[0/90]6, and quasi-three- dimensional Q3D [0/9O]6 weave specimens investigated experimentally in [13] are analyzed following the flowchart outlined in Figure 5.1. As can be seen from the damage sequence of the tested specimen shown in Figure 5.16, the damage is concentrated to a local region of about one unit cell of plain weave and quasi-three-dimensional weave at the center of the test specimen. Therefore, the shell-solid coupling technique can be employed for computational efficiency. 134 Deflection Figure 5.16 Damage sequence in Q3D samples from experiment [1 3] Rigid impactor 9.53mm in diameter Clamped along the \ edges Shell solid coupling 28.14mm 9.96mm/'\ ‘1 Symmetries along t e Local solid symmetric planes model Global elastic region modeled with shell elements S4R Figure 5.17 The generic shell-solid model used for L[0/ 90]6 , 2D[O/ 90]6, and Q3D[0/ 90]6 specimens 135 To model the laminated L [O/90]6 , plain weave 2D [O/90]6 , and quasi-three- dimensional Q3D [0/90]6 weave specimens, a generic shell-solid model as shown in Figure 5.17 is created. The test specimens are of 76.2mm x 76.2 mm x 2.04mm (length x width x thickness). Using symmetry, only a quarter of the plate of 38.1mm x 38.1 mm needs to be modeled. At the center of the specimen, which is 9.96m x 9.96 mm, a solid model is used while a shell model is used for the global region. The material properties used for the local solid and cohesive sections given in Table 5.2 are used here. The global shell model is the same for the three specimens including L[0/90]6, 2D[0/90]6, and Q3D[0/ 90]6 , except that different homogenized material properties must be used for the shell elements S4R to account for the three different microstructures. Table 5.4 shows the homogenized material properties computed using tow material properties of Table 2.6 while Table 5 .5 shows the homogenized material properties computed using the tow material properties given in Table 2.8. As explained in Section 2.7.2, the tow material properties in Table 2.8 are interpolated by scaling to ensure that the computed effective stiffness is the same as the measured stiffness from experiment of cross ply laminates. In Figure 5.17, the impactor is modeled as a rigid surface associated with a point-wise mass equal to 12.86 kg as used in the experiments. It has a spherical surface with a diameter of 9.53mm, which is the same as the impactor used in the experiments. An initial velocity of 3.1 m/s in the vertical direction is prescribed to the impactor, simulating the impacting velocity measured during the test. Additionally, a gravitational force in the vertical direction is applied to the impactor. 136 Table 5.4 Homogenized material properties using tow material Table 2.6 Material Err = 522 E 33 012 613 = 023 V12 V13 = V23 (GPa) (GPa) (GPa) (GPa) Laminate 23.58 9.08 4.14 3.75 0.09 0.36 [0/90]6 Plain Weave 21.25 8.64 3.81 3.54 0.10 0.37 [0/9016 030 Weave 20.73 8.64 3.77 3.66 0.10 0.38 [0/90], Table 5.5 Homogenized material properties using tow material Table 2.8 Material E11 = 522 E33 012 013 = 023 V12 V13 = V23 (GPa) (GPa) (GPa) (GPa) Laminate 13.8 4.84 2.35 2.08 0.09 0.36 [0/ 9016 Plain Weave 12.97 5.22 2.32 2.14 0.10 0.37 [0/90]6 Q3D Weave 12.69 5.20 2.31 2.17 0.10 0.38 [0/9016 To simulate the specimens used in experiments, the local solid models are modeled with continuum elements and cohesive elements as illustrated in Figure 5.18. 137 Interfaces modeled with Plies modeled with solid cohesive elements COH3D8 elements C3D8R (a) Laminated [0/ 9O]6 (only 2 cohesive layers are shown) Interfaces modeled with cohesive Tows modeled with . 4.65mm 0.66mm solrd elements (b) Plain weave [0/ 90]6' (analytical model at the left, illustrative model at the right) Interfaces modeled with Tows modeled with cohesive elements solid elements 0.66mm 4.65mm (c) Quasi-three-dimensional weave [O/ 90]6 (analytical model at the left, illustrative model at the right) Figure 5.18 The local solid models 138 In Figure 5.18, each composite ply is modeled with one solid C3D8R element through-the-thickness. For the irregular regions in the plain weave and quasi-three- dimensional weave structures, such as the undulation zones and matrix pockets, six-node linear triangular prism elements C3D6 and four-node linear tetrahedron elements C3D4 of reduced integration with hourglass control are used. Since the reduced-integration element formulation considers only the linear part of the incremental displacement field, the major drawback of this formulation is mesh instability known as hourglass modes. To minimize mesh distortion due to hourglass modes, the default enhanced hourglass control in ABAQUS is used for the solid elements. In ABAQUS/Explicit, localized stiffness reduction associated with damage can cause excessive element distortion, leading to numerical convergence difficulties and reducing the stable time step. These two factors may cause the analysis to run slowly or even to abort. In order to prevent excessive mesh distortion and allow the analyses to complete successfully, the elements are removed from the mesh once their dwage variables reach 0.99, which was widely used in studies of composites such as [9]. In this study, it was also found that element removal upon damage variables above 0.99 led to analysis termination due to excessive element distortion. Element-based surfaces are defined on the interior of solid bodies for use in modeling erosion due to element failure. The general contact algorithm removes contact faces and contact edges from the contact domain and, if an interior surface is defrned, activates newly exposed surface faces as elements fail. Thus, element-based surfaces are used to describe eroding plies/tows in this case. 139 The resin-rich interfaces between plies are discretized with COH3D8 and cohesive elements that model delamination by means of the traction-separation laws. To mesh the resin-rich interfaces in the undulation zones where COH3D8 elements can not be used, six-node three-dimensional cohesive elements COH3D6 are used In the study, the thickness of the cohesive layer within each layer is of 0.1% of the nominal ply thickness for the cohesive layers with undulation in Q3D weave and plain weave. For the interfaces between layers in laminated composites, cohesive layers with zero-thickness are used. The cohesive elements are removed to avoid uncontrolled distortion when the stiffness reduction reaches 0.99 [9], and they are replaced by the penalty contact formulation between plies specified in the general contact algorithm of ABAQUS/Explicit. Two algorithms for modeling contact interactions are available in ABAQUS/Explicit. The general contact algorithm allows the contact with eroding bodies and fiiction between the plies when the cohesive elements failed, while the “contact pairs” algorithm does not allow these features. Therefore, the general contact algorithm is used in this study. A mesh size of 0.25m x 0.25mm is used for the solid elements and the cohesive elements, and a mesh size of 1mm x 1mm is used for the shell section. Each ply of 0.17mm is modeled with one element in thickness direction. Due to the small sizes of the elements, the typical stable step increment calculated in ABAQUS/Explicit is found to be 10’8 second. Variable mass scaling is applied to the elements with incremental time step less than 10—8 second. Overall, the increase in the total mass of the model is kept under 1%. Fixed boundary conditions are applied along the edges of the shell, and symmetries are imposed to the associated symmetric boundaries of the quarter plate. The shell-to- solid coupling constraint in ABAQUS is a surface-based technique for coupling shell 140 elements to solid elements, and shell-to-solid constraints couple the degree of freedom such as displacement and rotation of each shell node to the average degree of freedom of the solid surface in the vicinity of the shell node. The shell-solid coupling constraints are used here to connect the global shell elements to the local solid model. 5.5.2 Result discussion The three models associated with laminate L[0/90]6, plain weave 2D [0/90]6, and Q3D[0/90]6 weave specimens are analyzed by ABAQUS/Explicit. The load-deflection curve from experiments [13] is given in Figure 5.19, while the load-deflection curves from simulation are shown in Figure 5 .20. All the load-deflection curves show that the Q3D woven composite can absorb more energy than its traditional 2D and laminated counterparts. As can be seen in Figure 5.20, the Q3D woven design obtained the largest peak load and the maximum energy absorption, while the laminated design has the smallest peak load and the minimum energy absorption. By comparing Figure 5.19 with Figure 5.20, it can be seen that the computed load-deflection curves show much larger energy absorption differences among the three different structures than those observed in experiments. While very close peak load correlation between experiment and simulation is obtained for laminated deign, the peak load difference for the Q3D woven design is as large as 50%. Furthermore, the simulation curves also show larger stiffiress than the experimental results. The discrepancies in the peak load and stiffness may be caused by the lack of the strain rate effect and plastic deformation in the material model employed because good correlation in the peak load and stiffness has been achieved in the static analysis of Section 5.4.3 using the same material model and solid-shell coupling technique. In the current dynamic analysis of high impact energy of 61.79 J (impactor 141 with a mass of 12.86 kg and a velocity of 3.1 rn/s), the strain rate effect may not be neglected anymore. This may also explain the small change obtained in the displacement curves in Figure 5.20 using different material properties. The computed maximum deflections for all the three designs are consistently smaller than the results from experiments. The load-deflection curve (Figure 5.20 (b)) using the less stiff material properties from Table 5.5 approaches closer to the load-deflection curve from the experiments than the one shown in Figure 5.20 (a), which is computed using homogenized material properties from Table 5.4 computed with the tow material properties directly form the supplier as given in Table 2.6. The strength and fracture toughness material properties used in the simulations are the same as shown in Table 5.2. The maximum deflection difference may be due to the following reasons: 1. Different material properties are assumed in computation than those used in the experiments due to the lack of the tested material properties. 2. The CDM model assumes no plastic deformation for the composites. This assumption may prevent the large deflections obtained from experiments. The largest deflection using similar assumptions was 6 mm for 4.4mm thick laminated plates subjected to impact loading [9], compared with 4 mm deflection for 2mm thick structures in the current study. 3. In this study, the stress based quadratic failure criterion (Eqns. (4.61) and (4.62)) is used for the prediction of delamination initiation under mixed-mode conditions. In the quadratic failure criterion, compressive inter-laminar normal stresses are assumed to not affect damage initiation. This may lead to the inaccurate, large delamination in the top 0/90 interface of the [0/90]s plate; and accordingly, the computed energy absorptions for 142 the three structures are consistently smaller than the experimental results. Results may be improved by replacing the quadratic delamination initiation criterion by the criterion proposed by Hou et al. [14], which takes into account the constraining effect of the compressive inter-laminar normal stresses in delaying the delamination initiation Figure 5.21 shows the comparison of the damage. As can be seen, the damage at the center of the specimens due to the penetration of the impactor is captured with one unit cell modeled with damage. Load (kN) N deflection (mm) Figure 5.19 The load-deflection curve from experiment [13] 143 P Load (kN) w P 9' o . Load (RN) 5” ' v I ' . 4:6 5.?) 7 37:0 Deflection (mm) (a) The load-deflection curve using shell material properties in Table 5.4 I v v — 20 -— Q3D ‘ ifo ‘ 4fo Deflection (mm) (b) The load-deflection curve using shell material properties in Table 5.5 220 Figure 5.20 Load-deflection curves from computation 144 ...‘vo‘nl-I'rcl ...»ucuv pa..~.‘-c..- - I .1.- .....nnq..v‘:\ :- < .- .... ....ouI-.o~o.~ .- -:-.. u-ru I....--.-..-n-:- - ...”..th \uo~. ... r 1.64%..“ r. 7 o.- ....uu 'IsI‘Il. v. I'l—Il v...-tl- u-..----~--.IOIII|-oa~- v-QIOIiO '0 nor-n .9 a o - v a o - e V v s 1 y . ‘ o I . -I-~a---uuc-- a... ~.. ...... 3 curl-onncivel.t we Ill-...hnt‘lI It! ... 4 - 2D [0/90], (0) Q3D[0/90]6 (a) Laminate L[0/ 90]6 ain weave (b) Pl Figure 5.21 Damage comparison (experiment [13] on the left, simulation on the right) 145 5.6 Discussion of the delamination initiation criteria For all the delamination studies in this work, the cohesive element formulation implemented in ABAQUS based on the original work of Camanho and Davila [15] is employed. This formulation uses the stress-based quadratic failure criterion Eqns. (4.61) and (4.62) for delamination initiation prediction and the energy based failure criterion Eqn. (4.63) for delamination propagation prediction. As pointed out in the current studies of the laminated [O/90]s and the work of Ayrnerich et al. [16,17], this formulation predicts inaccurately large delamination at the top 0/90 interface of [0/90]s laminates. Confusingly, by implementing the cohesive formulation proposed by Camanho and Davila in LSDYNA Hu et al. [18, 19] did not find the inaccurately large delamination at the top 0/90 interface of [0/90]s laminates in their studies. Ayrnerich et al. claim that the large delamination is caused by the delamination initiation criterion used in this formulation because the stress-based quadratic failure criterion, Eqns. (4.61) and (4.62), does not include the compressive inter-laminar normal stresses in the decision of damage initiation. Furthermore, Ayrnerich et al. claim that correct delamination at the top 0/90 interface of [0/90]s laminates can be predicted by using the failure criterion proposed by Hou et al. [14] given by < > 2 2 2 31— + -T—2 + :3— =1 for 1120 (5.1a) 2,16' 2.26 T36 2 2 2 f 2 2 7—2 + 2- —81 =1 for — Q—HLszl<0, (5.1b) T20 1'36 2'16 8 7224-1 2 No Delamination for 11<— —8——3—— (5.1c) 146 In the above failure criterion, the constraining effect of the compressive inter-laminar normal stresses in delaying delamination initiation is assumed. However, Ayrnerich et al. did not provide any information concerning the mix-mode formulation and calculation of the damage variables for the criterion proposed by Hou et al. [14]. The mix-mode formulation and calculation of the damage variables are critical to the implementation of this criterion because the damage variables control the initiation and propagation of the damage. In Camanho and Davila’s cohesive formulation, a mixed-mode relative displacement associated with the stress-based quadratic failure criterion, Eqns. (4.61) and (4.62), is defined as follows 6,, = «512) + 622 + 632 (5.2) where 6,-(1' = 1,2,3) is the relative displacement associated with Mode 1, II, and 111. To investigate the effectiveness of the delamination initiation criterion proposed by Hou et al., a mixed-mode relative displacement associated with the criterion by Hou et al. (Eqn (5.1b)) is proposed as follows 5,, = Jazz + 532 — 8<612> (5.3) Impact tests are performed on a two-layer [0/90] laminate model to study the delamination effect using the proposed formula. Only a quarter of the two-layer [0/90] laminate is modeled as shown in Figure 5 .22 (a) due to symmetry and it is impacted by an impactor of diameter of 9.53mm and mass of 2kg with an initial velocity of 0.5 m/s. To monitor the delamination initiation and propagation, damage is only allowed to occur to the nine cohesive elements at the center of the plate as illustrated in Figure 5.22 (b) and 147 the rest of the model is assumed to be elastic using the material properties from Table 5.2. The delamination initiation using Eqn (5.3) is found to be slower but the delamination propagation is much faster than using Eqn (5.2). Figure 5.23 shows the post-impact damage of the nine cohesive elements. In Figure 5.23, a value of zero indicates no damage to the element. The element damage increases with the value increases and a value of one represents failed elements. This study indicates that the inaccurate delamination predication problem may not be solved by only replacing the quadratic delamination initiation criterion with the criterion proposed by Hou et al. Modifications of the delamination propagation criterion in the formulation of Camanho and Davila should he considered. This is consistent with the conclusions of previous studies in the literature [20 21 22], in which it was found that the delamination propagation criterion is the dominant factor in delamination prediction. (b) Figure 5.22 A quarter of [0/90] laminated plate (a) with damage specified to occur only to the nine elements (b) at the center 148 Figure 5.23 Post-impact damage of the nine elements (a) using Eqn (5.2) (b) using Eqn (5.3) (0 denotes no damage; damage increases with the value increases) 5.7 Final comments of three-dimensional dynamic damage analysis The developed intra-laminar damage model combined with the cohesive zone model in ABAQUS/Explicit and the shell-solid coupling technique are used to investigate the performances of transverse impact resistance in laminate, plain weave, and Q3D weave structures. The results demonstrate the advantage of the Q3D weave over its plain weave and laminate counterparts in damage resistance. The computed maximum deflections for all the three designs are consistently smaller than the results from experiments. The maximum deflection difference may be due to the lack of plastic deformation assumed in the material model, the delamination initiation criterion used, and the material property difference between computation and experiments. The preliminary studies of the delamination initiation criteria show that the inaccurate delamination predication problem cannot be solved by only replacing the quadratic delamination initiation criterion with the criterion proposed by Hou et al. and thus 149 modifications of the delamination propagation criterion in the formulation of Camanho and Davila should be investigated. 150 10. ll. 12. 13. 14. REFERENCES . L. lannucci, M.L. Willows, An energy based damage mechanics approach to modelling impact onto woven composite materials—Part 11. Experimental and numerical results, Composites: Composites: Part A 38 (2007) 540—554. L. lannucci, M.L. Willows, An energy based damage mechanics approach to modelling impact onto woven composite materials—Part 1: Numerical models, Composites: Part A 37 (2006) 2041—2056. L. lannucci, R. Dechaene, M. Willows,and J. Degrieck, A failure model for the analysis of thin woven glass composite structures under impact loadings, Computers and Structures 79 (2001) 785—799. R. Krueger, J. G. Ratcliffe, and P. J. Minguet, Panel stiffener debonding analysis using a shell/3D modeling technique, Composites Science and Technology 69 (2009) 2352—2362. R. Krueger, J. G. Ratcliffe, and P. J. Minguet, Analysis of composite panel-stiffener debonding using a shell/3D modeling technique, NIA Report No. 2007—07, NASA/CR-2007-214879; 2007. . R. Krueger and P. J. Minguet, Analysis of Composite Panel-Stiffener Debonding Using a Shell/3D Modeling Technique, NASA/CR-2006-214299, NIA Report No. 2006-02. R. Krueger and P. J. Minguet, Analysis of Composite Skin-Stiffener Debond Specimens Using a Shell/3D Modeling Technique and Submodeling, NASA/CR- 2004-212684, NIA Report No. 2004-04. ABAQUS/Explicit User’s Manual, online version 6.9-EF1 (2009), Hibbit, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. C.S. Lopes, P.P. Camanho, Z. Giirdal, P. Maimi, and B.V. Gonzalez, Low-velocity impact damage on dispersed stacking sequence laminates. Part 11: Numerical simulations, Composites Science and Technology 69 (2009) 937-947. B. G. Gulker, Energy Separation in Impacted Composite Plates, M.S. thesis 2009, Michigan State University. J. Fish et al. Elastomer Impact Simulations, Jacob Fish et.al, online paper, http://www.scorec.rpi.edu/REPORTS/2005-1 1.pdf. S. Kamiya, H. Sekine and Y. Yagishita, Computational simulation of interlaminar crack extension in angle-ply laminates due to transverse loading. Journal of Composite Materials 32 (1998), 744—765. K. Rosario, Impact Properties of Quasi-3D Composites, M.S. thesis 2008, Michigan State University. J .P. Hou, N. Petrinic and C. Ruiz, A delamination criterion for laminated composites under low-velocity impact, Composites Science and Technology 61 (2001), pp. 2069—2074. 151 15. l6. 17. RP. Camanho and CG. Davila, Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials, NASA/TM-2002-211737. F. Aymerich, F. Dore, and P. Priolo, Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements, Composites Science and Technology 68 (2008) 23 83—2390. F. Ayrnerich, F. Dore, and P. Priolo, Simulation of multiple delaminations in impacted cross-ply laminates using a finite element model based on cohesive interface elements, Composites Science and Technology 69 (2009) 1699—1. 18. N. Hu, Y. Zemba, T. Okabe, C. Yan, H. Fukunaga, A.M. Elmarakbi, A new cohesive model for simulating delamination propagation in composite laminates under transverse loads, Mechanics of Materials 40 (2008) 920—935. 19. N. Hu, Y. Zemba, H. Fukunaga, H.H. Wang, and A.M. Elmarakbi, Stable numerical simulations of propagations of complex damages in composite structures under transverse loads, Composites Science and Technology 67 (2007), pp. 752—765. 20. A. Turon, C.G. Davila, P.P. Camanho, and J. Costa, An engineering solution for mesh 21. size effects in the simulation of delamination using cohesive zone models, Engineering Fracture Mechanics 74 (2007) 1665—1682. G. Alfano and M. A. Crisfield, Finite element interface models for the delaminationanalysis of laminated composites: mechanical and computational issues, Int. J. Numer. Meth. Engng 2001; 50:1701-1736. 22. J. Chen, E. Ravey, S. Hallett, M. Wisnomc, and M. Grassi, Prediction of delamination in braided composite T-piece specimens, Composites Science and Technology 69 (2009) 2363—2367. 152 Chapter 6 Conclusions and recommendations 6.1 Conclusions 1. Using the developed numerical homogenization method, good correlations with experiments are achieved for the comparison of stiffness between Q3D woven composites and their laminated and two-dimensional woven counterparts. Parametric studies using this homogenization method show that undulation and fiber volume fraction are the two key factors affecting the effective in-plane stiffness of woven composites. Based on the results of the parametric studies, guidelines are presented for designing Q3D woven composites with high in-plane stiffness. Furthermore, the developed method can be used to perform homogenization of composites with complex geometries. 2. Using the cohesive zone model in ABAQUS/Standard, double cantilever beam (DCB) tests are performed on composites made of unidirectional layers and plain weaves. The effects of stiffness on delamination initiation and propagation in these composites are discussed. Both DCB and end notch flexure (ENF) tests are modeled for woven composites with cohesive layers inserted at different interfaces. The undulation effects on delamination propagation in woven designs are investigated. An important conclusion from the ENF test is that delamination resistance can be improved by eliminating the horizontal pure matrix interfaces in laminated composites and using interlocking layers as in the Q3D weave designs. 3. A three-dimensional continuum damage mechanics model for predicting the initiation and propagation of intra—laminar damage mechanisms is developed by integrating the stress-based failure criteria for damage initiation, fracture- 153 mechanics based criteria for damage propagation, and the crack band model for alleviating mesh sensitivity. Thermodynamics irreversibility in this model is ensured as a result of the positive evolution of damage variables due to consistency conditions. The three-dimensional damage model is implemented in ABAQUS/ Standard using a user-written material subroutine (UMAT). The effectiveness of this model is demonstrated based on the good correlations between numerical results and experimental results reported in the literature for laminated plates subjected to transverse static and low-velocity impact loadings. The progressive damage process of fiber-bridging in small-angle laminated composites subjected to transverse low-velocity impact is successfully simulated. The results demonstrate the necessity of in-ply matrix cracking for the successful simulation of delamination and fiber-bridging in laminated composites subjected to transverse loadings. . To study the impact resistance of composites, composite plates made of unidirectional layers, plain weaves and Q3D weaves are investigated in ABAQUS/Explicit. With the use of the shell-solid coupling technique, both numerical accuracy due to three-dimensional solid elements and computational efficiency owing to shell elements can be achieved. In the shell-solid model, material properties computed using the homogenization technique developed in Chapter 2 are assigned to the shell elements in the global region. The damage of local region is modeled using the developed three-dimensional continuum damage mechanics model. A user-written material subroutine (VUMAT) for intra-laminar damage is implemented in ABAQUS/Explicit. The results demonstrate the 154 advantage of the Q3D weave over its plain woven and laminated counterparts in damage resistance to transverse impact. 6.2 Recommendations 1. It is recommended that a mix-mode formulation be derived by using the criterion proposed by Hou et al. [1] for delamination initiation prediction and modifying the delamination propagation assumptions in the formulation of Camanho and Davila [2]. In the development of the continuum damage mechanics model, it is assumed that strain-softening is solely due to the degradation of material stiffness without any inelastic behavior. After the material is damaged, unloading leads to complete closure of the cracks. For fatigue analysis of composites, plastic deformation might be needed. The assumption of lack of plastic deformation might prevent the correct prediction of the large deflections obtained from experiments in the impact analysis of current study. The developed three-dimensional continuum damage mechanics model assumes no strain-rate effect in the material behavior. Good correlation with experiment is achieved for static tests. In order to obtain good correlation for impact tests, the strain rate effect should be considered. The modeling of woven composites is very complex when using continuum elements for intra-laminar damage and cohesive layers for inter-laminar damage. Delamination modeling may be integrated into the three-dimensional continuum damage mechanics model to ease modeling efforts. 155 REFERENCES 1. J .P. Hou, N. Petrinic and C. Ruiz, A delamination criterion for laminated composites under low-velocity impact, Composites Science and Technology 61 (2001), pp. 2069—2074. 2. PP. Camanho and CG. Davila, Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials, NASA/TM-2002-211737 156 I I I I I I I I I