v .191 999‘999‘9.C ‘ 99 :999}:V , 9, 9, .199 95\'|‘ “Yaqllfi W 5‘ 9’,9 «999 9, l 9. 9’4 . 4 5 ‘ . “ ' '. ‘ I 55 'I "—‘ ‘1 n' —'| ' 1 9.55%‘915‘9‘15 'f h 'b‘ ' 5". , 4. . I 9 v "I, , 9 99\y.9 ‘ ‘,"95 95 'i' .4929 ,.;,.9 9' _ 5', 9 . . d .9,‘ 5:9'.1J “ V5! ”991995 ‘9 :59 ,1,999'5~ ‘ 99:94.," “99'2".99,’ 9’99 9 959 99 1:. ,L, "9“; 99 ‘99‘9'59392'1'9-9 1'1 , , $59: '9955 999 'ZQJ'U‘)“ 91,, 5599995 514.93.: 99 919995991, " 5 ,9‘9“ 1,,29'5 ~". “.5199. .9 1‘ "'1' ,99999' ‘ififi'kh 5‘ ‘ 9‘ 9 9 ' 599;" 93:19:59. 9," ‘9 93:13? ' 99.9” 9,-9.9, 5"“= 91:99:99 9195‘? ‘9‘9-9 “ W ‘ 9 ‘9. 93.9% E“. 9 EU' {9‘9 999554919999 555999; 4 99,9 91ft] 9'5 9999.9 599954951959 J“); “9999999999- 9:1?‘9'1‘999 , , 9. "'5‘91'i9 9:9 9: 99 99999 9 ,9}. 2’ 5 9‘91‘L19.‘99 ‘599I199959999, 9 3 919‘!,9‘{9..'{99,: 9,“..99939999 5999 991599 999999 9199 1‘39} . . 999 (9999;, ' ‘ M59W1‘M'1 9 fimhflfinawd .w.%.‘.flfim.9gg{ww5 99‘“ 9’59‘13995‘ 59‘19995‘1119599‘ ‘9"; "999 ‘9 '9‘" '9' ‘9‘ " 9,9,.‘91-99 = 5353‘ ' " . 1753‘? 93993 999.259 591199329: 9, . ,55"'!995 9'99‘9‘9999 959 ”W {{ 949115‘ 99‘, J .. . 19 999 9:479";‘.'.95{‘+:“9£h5 99-11? 35:93.95}4.9f143954af';9353:;9FIC9'1fi556‘ir9951fi‘i1"5‘3F-95g -, l 955‘ 1 1‘5 :9 9'95' ‘ 1 19'95‘1191""'1 '59‘51): .9: ,9“‘59 '9'1111-‘7.‘5""519‘"1'i9."'459 -1 1- 9 99 59. 9 999 59 "9999199 9999,.) 9 9 9 9 9 9'1"} (f. 99 99 f I“E':‘-.;*f‘—,‘-'Vf‘ t9: 9 _ _9 55.9195 51591991915 . 5N {9993.9999" 9::e 1"11995'1159599jni I595 5 -:9E5!‘{95{:9‘91,.I:I19’;‘9 21‘959555‘y3. 9939.19.49 99959:, 513111559919," '999 ff; ‘2:;:95_£99§9,9993§§.r9- 9.9. 15‘ .15555'55""11"113":t11'11'1""'11"¢ fli‘ “.551 ' *‘3‘9 $3915?“ (:3 ' £99.39 ' »’ 13 '9 '1‘,“ 9999 5999992951.“.Q1995 95"-9',5:'9)1;9‘ ,E E999;9 59 9mm. 9""91‘5‘9 9,9!;'91999¥}4995515.§:: ,‘9 9525:} 95 9596,!) 5;;‘993959 15355429955} ..9 9.9 , .2 ‘ ‘ 9 . . .4, V t: I, ‘.4 ‘!¢ .. 3'. 9 ‘. '.‘ .‘ 4999,1594} 5'9919919l1951' 59199, 99.9‘, , ,- ‘19 9: 7:5‘9.35{'19.‘;‘_.,?9:9:2‘Ii: 99.99:". , M',“'9fi,:‘.~;5."3j.'9.'iq:h!,$739, ,"z "9;!" J; 99‘199 5999.15 5999 9 ‘i55 9.99999 9 ’ 991999599 9 999‘, £53553, , 9919999; 99,993,999 999— 9’, {1'9 39 :L. .29 ‘ 5.13:,9'9" "91 ' '99 "5. {9| ,9, 99;.- s ' .9 ‘. 5W5 " 9 )9“ “9‘?“ "- .9 m9; » ‘ 9' 1“ 9:99 9 u I ' 5" O '5 5‘) 59: ‘ "996‘ , . 9. 9;:(E‘9‘95‘9' fl99551‘91919992. ‘555' :5%’1'%'299.99939991 9999.999 99.999.91.999; 9‘ .9 59) 9'39995'1' 9"5'29’9 ' ,, ‘ '1.':‘I999 9“"59.‘- v:9 9x9 9. _ .9 5" 99.5. '5‘". I 19“.! 5 :1 1‘ ‘5 ‘. ‘ '9'”; .6541. 99 4.";2 ‘ " 1,.1'71‘714 ,o 9 591i "9‘13; ~9v’5 ‘ 'V' 1.! (9L "1'1""111111"EW'9"'1'""1:"".1"1‘1""'1" 5 l I“ M‘ "11' ‘1 ‘ 11"; '3 .591W‘94' fl: ”4'52;";1';g 9,9; ' 9."'.’"" 5‘53}: .‘i';'9 .. 9 9 " :‘,99?;.'.‘;“"§9;95 9231‘ 9- .99.; 4. 9.399 9 9Wm~“”‘v',9,u .. ,9.. ' 99!: 993991,, 9': 992‘: 9 9999 -' 22‘: .999 9.9:- 3.; 5‘LL1'H'E191'91'1‘5W‘ 91" '1'" '1 “1‘ 999-. 1' ‘519'95 "5 9 95411.5‘56‘1"““'*'55i. ;;§.52-55595;11'9 ' '5 I 9515?“? 2'99‘91;11'" , H.151”; "9" 9 9,, 3."" ‘gt -1‘ 1|_.1:\- ~‘j'9'2193'9'59: 9115'“ 193%“? . I “n,¢.~wmau.u9;mmvwww' 9 -‘9 w gw'~9~~;9.-- ,5 ‘5“.4" .191'9515‘.l"|97 '9 5 -11 '1 ‘ ‘ 3,99,.99 45595: 14.5.1; 91.294511. "'95‘9:'i. 9 ,"J‘: .. ' ' 29; 3 9,933.. T192511 ‘ 19, n9 1'. u .‘ nay ' 29919 1 '1" If 99;' ;5 3.; i5; 9" 5’ 45. .95 '9. , 5" .1 , 2; 9': ‘ .' '5' , , I 4 9 7|: 9' 9 .591. n ' ' 99 9' 9 f 5;!) 99 I .'91'1l;:9'551 9 9" In); 5'9‘91‘5'”: 94' v.5 ' ,9?“ . 9 9 ‘ 9 9- 9 599;, '9‘ ' 4:49.99‘1'H" ‘ 1 . 5 , 1.51' ;1'.9.".9L $1.959 145.9195 '99 1:“ '°""':"‘.14'9 ‘ '9‘5 5‘51.) 9." 9' . '.-‘.9" ' '99:'51." '9"9.1'1f"9‘4'.49¢ .. 99 9 9 .9, - - 9. I v‘ "t. 9 9:'-‘-‘-". 91?! "t ' 5395-“ ' 2".99'1 " , “99' 9 .;.9-‘1:99 ‘h‘ 9.39;" 9‘9“. "9.??‘3 ~'9, {95:19 .‘ 1 99 :, “112 ‘9' 991'""‘," " ‘ 9 9 . ‘ ' ,9‘9 I" 99 9999" . 9-i; 93" -‘ '1 91.1"".15“ ‘22": .9199 “9.99999 9 999,9 9.91.99.99, 99,9 '9 .I , , 9 , 9.1- . 5| ,‘993‘99995 9 '99 {1, 9 59'913'59‘7'93'5‘99'W'1N” 0999559; "'1 '91‘ ". |"5951'1‘¢i9"9)‘9:1,‘5VI";.1u1t-“9'9‘9'59I ‘ 9 ':;}1\555'91'9995"59 9-9 M9 9 9 m9 9 9|“ 5 9"‘9.: 9-, “9'99. " 5‘9 ',' .55'59'1‘ 1 1‘1" '1'11 9.99.5.1, '1“ 1 '19 .‘I' I' ‘9 9' '4.9,9"99 99,9'9199 942'”. 9. 9999'6" . I ‘ 999‘9 19,‘9'."1.9'; .9999. 9,9 ‘9 ,9'. 99 9, 99999 99.99.959-99-9 ‘99 9 '1‘999’9 29;)159912559f29'19195955999 9 49951.1"'5"11"'.:5."09"-115:'9.:'11"‘§9.:"159‘ A 09:5‘) 1"1'911’Ho 9 99! :Ei'f1z1997.'9109'9;9"99 '35) 5.5;“. , 'P95‘50‘, 9':1‘,,.99 ‘ 99: ":1. 5"" ' "' "'55'9‘1'1 “01"1‘5' ll11"95l ' 9 "i ' . H" ‘5'! ‘ [95"‘5: “5 II 5 9153551951!“ 9 995. .49 ' {,59m9‘fi.”9&999‘9959'9fl‘www99‘ 999959999999‘9mw .Aflw,flfl%5;§ “' ‘" 9‘99 “"3 '1 ‘99 ... - 9‘9 !9i; 9:.” ‘ ‘ "9‘52? 3' 999.“:9?!‘ 99 9:999"9 99.. '-‘.,' .9 '9911‘5}!|5;5’99'lf559‘1‘55;5,115I99999999“9'9“995,99999,!93995‘999,999957,955,9599 “435959395959 5‘5599 9‘99:.~,99...159;;9 99995;; 95555559 WWI; 1‘ 5559"“ 55 9 99 99. 9 ”'9, 9 5,9 599 ~45 [99,9 , 99) 9'99 9 5"1'" 3,} .9 ~ "“"‘)~, 9999 9999499 9 5.:,{ 9599319. 9 99.5 ,39 999 9, ,. , 9 ‘99 9 99.9- 19 9 i 999 9 9 9i ‘9 ‘5999‘9‘9‘91 t? = '992‘9999 .1. ‘91.... '9': “9{'1“"""59"§‘!{ELI-9' l 9.5 . . 595351.“ 115"? 91949999999959: -, 999599, 9,, 99.999 99,999 99999 99999199 9999; 9-99 ~ 9499 “5'52‘13431‘1-99fi125‘4599'E-‘W " ”l’ " 999 ,999‘9‘9 ‘99 9999 - . 9.99, ,9 . 99,. . , ~ 9 .:91)55"51,9 ,, 5 ,9 ' 9 95 5519999§,91?9"599 39:15:95 9 '9, 95 ’ 9. 9595 ‘1 2’9 '9‘ . 9 9 1 . ‘5' ,1 5‘ 55‘ 9 1r; . 5,9.5' i5 9 9199 {9 9 9999.99 9‘ 999 9 95‘] ,5 )9 79‘ 959.995. 97.159539”? " 599999 99 9.9.99“; 9‘9 GAN STATE IVERSITY LIBRARIES lllllulll llllllm | l l 3 1293 00885 0939 This is to certify that the dissertation entitled ANALYTICAL AND NUMERICAL STUDY OF PHASE TRANSFORMATION WITHIN A BAR COMPOSED OF A STRAIN-SOFTENING NONLINEAR ELASTIC MATERIAL presented by Jiehliang Lin has been accepted towards fulfillment of the requirements for Ph . D . degree in Mechanics mafia? M jO professor Date Z/Z/93 MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY University Michigan State PLACE ll RETURN BOX to move this checkout from your record. TO AVOID FINES return on or baton duo duo. DATE DUE DATE DUE m DATE DUE __Jl __JL___ :1? I ll MSU Is An Affirmative Action/Equal Opportunity Institmion “ottoman-9.1 *7 m _____fi_ ANALYTICAL AND NUMERICAL STUDY OF PHASE TRANSFORMATION WITHIN A BAR COMPOSED OF A ST RAIN-SOFTENING NONLINEAR ELASTIC MATERIAL By Jiehliang Lin A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Materials Science and Mechanics 1993 ABSTRACT ANALYTICAL AND NUMERICAL STUDY OF PHASE TRANSFORMATION WITHIN A BAR COMPOSED OF A ST RAIN-SOFTENING NONLINEAR ELASTIC MATERIAL By Jiehliang Lin This study is concerned with a model describing phase transformations within a strain-softening bar induced by stress. In an analytical study we consider two materials with particular strain-softening constitutive laws. These materials have a non-convex strain energy density function and their equilibrium equations may be non-elliptic. The parts of the bar with strains corresponding to the strain-softening part of the constitutive law are regarded as being in a different phase from the rest of the bar. Two loading de- vices are the main concern: one is the hard device in which the bar is firmly fixed at one end while the other end is subjected to a given displacement; the other is the soft device in which the bar is also firmly fixed at one end while the other end is subjected to a given dead load. A detailed characterization is given for the associated nonunique solutions to this problem. Here these nonunique solutions represent different ways to distribute phase transformations within the system. Energy functions for the solutions are also given. For both device loadings the problem is also studied within a quasi-static setting (which neglects inertial effects). Two separate criteria to resolve the nonuniqueness of the solutions to the problem are considered: a kinetic relation criterion and a minimum energy criterion. In the numerical study we consider the finite element method. It is shown that the nonuniqueness of solutions to the problem in the continuous system is translated into nonuniqueness of solutions in the numerical system. The FEM is then successfully im- plemented for each of the two separate criteria mentioned above. Copyright by Jiehliang Lin 1993 and my m y ACKNOWLEDGEMENTS It has been the author’s privilege to have had the following gentlemen serve on his doctoral committee: Professor Nicholas J. Altiero Professor Phillip Duxbury Professor Dahsin Liu Professor Thomas J. Pence From the beginning of my graduate education to the conclusion of my research and from conception to completion of the work for this dissertation has required not only my own devotion but also that of an experienced and patient advisor. I wish to express my sincere appreciation to Dr. Pence for his guidance and support during my research. I wish to thank my parents for their constant support and encouragement. Special recognition is due to my daughter, Michelle, for her tolerance of the demands a doctoral program places on a father’s time and disposition. And to my wife, Lingyuh Chem, where words are truly inadequate to express my appreciation for her continual encour- agement and understanding. I acknowledge financial support by US. Army Research Office under contract DAAL03-89-G-0089, and by the Composite Materials and Structures Center, Michigan State University, under the REF program. Table of Contents 1. Introduction ...................................................................................................................... 1 1.1 Statement of problem ................................................................................................. l 1.2 Literature review ........................................................................................................ 2 1.3 Section overview ........................................................................................................ 8 Part I: The continuum mechanical model for phase transformations in a strain-softening bar 2. The material model ........................................................................................................ 13 3. Two loading devices ...................... 17 3.1 Soft device (Figure 3) ............................................................................................... 17 3.2 Hard device (Figure 3) ............................................................................................. 18 4. Two phase equilibrium configurations with a single phase boundary ........................... 20 5. The solution region with a single phase boundary ........................................................ 22 5.1 The boundary of (2,2 ................................................................................................ 25 5.2 Lines of constant 10 in (2,2 ....................................................................................... 26 5.3 Lines of constant 0 in 012 ........................................................................................ 27 5.4 Curves of constants in on ...................................................................................... 27 6. The hard device energy function on $212 ........................................................................ 35 6.1 The hard device energy function on the boundary of £212 ........................................ 35 6.2 The hard device energy function for (212 on the lines of constant a ........................ 38 6.3 The hard device energy function for (212 on the lines of constant 70 ....................... 43 6.4 The hard device energy function on the curves of constant s .................................. 46 6.5 The minimum energy criterion for the hard device .................................................. 48 7. The soft device energy function on (2,2 ......................................................................... 51 7.1 The soft device energy function on the boundary of 9,2 ......................................... 51 7.2 The soft device energy function for (2,2 on the lines of constant c ......................... 53 7.3 The soft device energy function for £212 on the lines of constant 7,, ......................... 54 7.4 The soft device energy function for £212 on the curves of constant 5 ....................... 57 7.5 The minimum energy criterion for the soft device ................................................... 58 8. Quasi-static motions with a single phase boundary ....................................................... 59 8.1 Definition ................................................................................................................. 59 8.2 Dissipation and driving traction ............................................................................... 62 8.3 Admissible direction of movement of the phase boundary ...................................... 67 8.4 Comments on the admissible paths .......................................................................... 69 9. The kinetic criterion ....................................................................................................... 76 9.1 The governing differential equation for phase boundary motion ............................. 77 9.2 Phase boundary kinematics ...................................................................................... 79 9.3 Rate dependence of the admissible paths governed by the kinetic criterion ............ 81 Part II: Finite element approximation of the continuum mechanical model 10. Finite element analysis ................................................................................................ 85 10.1 Variational formulation of the equation over a typical element ............................ 85 10.2 Interpolation functions and element equation ....................................................... 87 10.3 Assembly of element equations ............................................................................. 89 10.4 Imposition of boundary conditions ....................................................................... 90 10.5 Iterative solutions to the system equation ............................................................. 95 11. Solutions of the system equation with N=2 ................................................................ 98 11.1 Formulation ........................................................................................................... 98 15_ N1 15.; 15,; 1: 11.2 Direct search for FEM solutions ........................................................................... 99 11.3 The hard device energy function A as a function of the possible nodal dis- placements .......................................................................................................... 106 11.4 Iteration solutions and the initial guess displacements ....................................... 107 12. Solutions of the system equation with N=3 .............................................................. 110 12.1 Formulation ........ 110 12.2 The hard device energy function A as a function of the possible nodal dis- placements .......................................................................................................... 1 11 12.3 Iterative solutions and the initial guess displacements ........................................ 118 12.3.1 Material A ..................................................................................................... 119 12.3.2 Material B ..................................................................................................... 123 13. FEM solutions for N>3 ............................................................................................. 124 14. Numerical results using the finite element method for material A ........................... 129 14.1 Minimum energy criterion .................................................................................. 130 14.2 Kinetic criterion ................................................................................................... 132 14.2.1 Several algorithms to find s (tM) .................................................................. 133 14.2.2 Computation of an “exact” solution path (o(t), 70(0) which satisfies the kinetic criterion ............................................................................................. 136 14.2.3 Solutions obeying the kinetic relation obtained by using the finite ele- ment method ................................................................................................. 138 14.2.4 Three ramp loadings with different loading rates ......................................... 140 15. Numerical results using the finite element method for material B ........................... 147 15.1 Minimum energy criterion .................................................................................. 147 15.2 Kinetic criterion ................................................................................................... 147 15.2.1 Computation of an “exact” solution path (0(t), 70(0) which satisfies the kinetic criterion ............................................................................................. 148 15.2.2 Solutions obeying the kinetic relation obtained by using the finite ele- vii ment method ................................................................................................. 149 15.2.3 Three ramp loadings with different loading rates ......................................... 152 16. Conclusions and recommendations for future work ................................................. 155 16.1 Conclusions ......................................................................................................... 155 16.2 Recommendations for future work .................................................................... ~..155 Appendix-A: Material A ................................................................................................. 157 A.1 The solution region (212 for material A ................................................................ 159 A.2 The hard device energy function A12 for material A ........................................... 160 Appendix-B: Material B ................................................................................................. 163 B] The solution region (212 for material B ................................................................ 165 B.2 The hard device energy function A12 for material B ............................................ 167 Appendix-C: Ordering of the solution paths in Figure 33 .............................................. 169 Appendix-D: Local analysis of the solution paths in Figure 33 near t=0 for material A. ............................................................................................................. 171 List of References ........................................................................................................... 174 viii Fi Fit Fig Fig Fig Fl g1 Fl gt FigL Figi FigL Figu - List of Figures Figure l. The stress response function of the bar with positive strain. ............................. 13 Figure 2. The non-dimensionalized stress response function of the bar which is used in this thesis. ................................................................................................... 14 Figure 3. Undeformed and deformed configurations of soft device and hard de- vice. ................................................................................................................. 17 Figure 4. The solution region f)“. .................................................................................... 22 Figure 5. The solution region 922. .................................................................................... 23 Figure 6. The solution region (212. .................................................................................... 23 Figure 7. The solution region 92,. .................................................................................... 24 Figure 8. The solution regions £21,, S212, 021 and (In. ...................................................... 24 Figure 9. The curves of constant s=(0.05, 0.1, 0.15,..., 0.95) in the solution region 912 with the choice of the function 3,,(7) given by (a. l) of Appendix-A. The limiting case of 5:0 corresponds to (ho-”(70) and the limiting case of s=l corresponds to both o=yo with (Doc and 0:6,, ................................... 28 Figure 10. It is con'vient to reorient the solution region of Figure 4-Figure 9 so that the constant s curves display a minima. We shall use this orientation for the solution region 5212 in the rest of this thesis. ............................................. 29 Figure 11. The $2 curve in the solution region 5212 (again the choice of the function 6,,(7) is given by (3.1) of Appendix-A). .......... . .............................................. 33 Figure 12. The level curves of constant hard device energy function A12 in the so- lution region (212 with the choice of the function 3,,(7) given by (a. l) of Appendix-A. ................................................................................................... 36 Fi; Fig Fig Fie Fig Fig FigU Figu Figure 13. The function cp,(o) with the choice of the function 3,,(7) given by (a. l) of Appendix-A. Note that (63)1 ensures that this curve is above the dot- ted line. ........................................................................................................... 40 Figure 14. The first derivative of function op,(o) with the choice of the function 6,,(7) given by (a.1) of Appendix-A. .............................................................. 43 Figure 15. The hard device energy function A12 in the solution region 5212. .................... 47 Figure 16. The minimum energy solution for the hard device. ........................................ 49 Figure 17. The loading-deformation diagram for the hard device under a minimum energy solution criterion. ................................................................................ 50 Figure 18. The level curves of constant soft device energy function 812 in the solu- tion region (212 with the choice of the function 5,,(7) by (a. 1) of Appen- dix-A. .............................................................................................................. 52 Figure 19. The minimum soft device energy solution on the constant 70 lines for 70 > 1. .................................................................................................................. 57 Figure 20. The minimum energy solution path for the soft device. .................................. 58 Figure 21. The loading diagram for the soft device under the minimum energy cri- terion. .............................................................................................................. 58 Figure 22. Some possible quasi-static motions in the solution region 012 as 70 in- creases from 70:0. ........................................................................................... 61 Figure 23. The functionsfnw) with the choice of the function 3,,(1) given by (a. 1) of Appendix-A for material A and by (b. l) of Appendix-B for material B. ..................................................................................................................... 66 Figure 24. The first derivative of the functions 0312(0) for the example function 6,,(7) given by (3.1) of Appendix-A for material A and by (b.1) of Ap- pendix-B for material B. ................................................................................. 66 Figure 25. The admissible directions for the quasi-static motions in the solution re- gion 912. The constant 3 curves are for the function 5,,(7) given by (a. 1) Fig Fig? Figt Figu Figu Figlm Figure of Appendix-A, but the diagram is representative for any function obey- ing (5). ............................................................................................................. 68 Figure 26. The admissible directions for the quasi-static motions in the solution re- gion 9'21- The constant 3 curves are for the function 3,,(7) given by (a. 1) of Appendix-A, but the diagram is representative for any function obey— ing (5). ............................................................................................................. 70 Figure 27. Admissible constant 0 paths that reach the boundary 70:1“,(0’) for the soft device. All points on the constant 0 lines can reach points on this bound- ary. .................................................................................................................. 71 Figure 28. Admissible constant 0 paths that reach the boundary 10=F,,(6) for both the hard and soft device. No points, other than the boundary points themselves, can reach the points on this boundary. ........................................ 72 Figure 29. Admissible constant 7,, paths that reach the boundary 70:1“,(0) for the hard device. Those points which are on the constant 7,, line with 5(y,)> 392%) can reach this boundary. ..................................................................... 72 Figure 30. Admissible constant 7,, paths that reach the boundary 0: CC for the hard device. Those points which can reach this boundary are on the constant 7,, line with oc< E's-(yak 6232(7) whenever oc< yo<1, or with oc< 3(70)< 3,,(70) whenever ya 2 1. ................................................................................... 74 Figure 31. Some admissible possible quasi-static motions in the solution region 912 as 7,, increases from 70:0. ............................................................................... 74 Figure 32. The hard device loading condition (149) for different values to. Here 0 = tom) < (00,)< < {0(6) < (0“) = 0° 80 that 0° = (10(3) > (100,) > ... > 00(6) > Figure 33. The relation between the solution paths corresponding to the loadings given in Figure 32 for the hard device. ........................................................... 83 Figure 34. The coefficients in Pascal’s triangle gives the number of phase permuta- Fig Figt Figu Figu Figm Figur Flgur FigUlt Fngrs tions at fixed phase volume fraction for two-phase problems with differ- ent number of elements. Here 2” gives the number of possible phase arrangements, and N+1 gives the number of distinct volume fractions of phase-I to phase-II. ......................................................................................... 94 Figure 35. The regions of four cases of possible 32 considering a given 70 for mate- Figure 37. The final trace of the different functions 66:2) for 70:08 shown in Figure 36. There is one solution to the composite trace 9(34'2)=0. ........................... 103 Figure 38. The different functions 9672) with two equal element lengths and 70:2 for material A. ............................................................................................... 103 Figure 39. The final trace of the functions 9(32) for 70:2 shown in Figure 38. There are three solutions to the composite trace 9(Hz)=0. ...................................... 104 . Figure 40. A schematic diagram of the possible solutions 52 for a given 70. .................. 105 Figure 41. The possible N=2 FEM solutions as represented in (212 for 70:2 and for 70:8 for material A. ...................................................................................... 105 Figure 42. The hard device energy functions Kay) with two equal element lengths and 70:08 for material A. ............................................................................. 107 Figure 43. The hard device energy functions X52) with two equal element lengths and 70:2 for material A. ................................................................................ 107 Figure 44. The hard device energy functions KG?) with two equal element lengths and different 7,, for material A. ..................................................................... 108 Figure 45. The hard device energy functions 1364}, 33) with three equal element lengths and 70:08 for material A. A similar picture results for material B at 70:08. ...... 112 Figure 46. The hard device energy function AG,» 53) with three equal element xii lengths and 70:2 for material A. Material B gives a similar picture. This picture is representative for all 70>1, both for material A and for material B. ..................................................................................... ' .............................. 113 Figure 47. The contour diagram of the hard device energy functions A52. 133) on the domain (219) with three equal element lengths and 70:08 for material A. A similar picture results for material B at 70:08. ................................... 113 Figure 48. The possible N=3 FEM solutions as represented in (212 for 70:2 and for 75.8 for material A. Comparing this diagram to Figure 41, it is noted that now two internal constants curves (s=1/3, s=2/3) are available for non-pure FEM solutions. .............................................................................. 115 Figure 49. The contour diagram of the hard device energy functions Adi-2, 'u'g) with three equal element lengths and 70:2 for material A. This picture is rep- resentative for all 70>1, both for material A and for material B. .................. 115 Figure 50. The contour diagram of the hard device energy functions NE, 33) with three equal element lengths and 750.97 for material B. This is different from the contour diagram for material A with 75097 (which is similar to the diagram in Figure 47). ........................................................................ 116 Figure 51. The possible N=3 FEM solutions to (217) as represented in 012 for 750.97 for material B. In addition to the pure phase solutions, FEM so- lutions are associated with the curves s=l/3 and s=2/3. There is no inter- section of 70:0.97 with 5:1/3 (dotted). There are 2 intersections of 75097 with s=2/3. each of which corresponds to 3 phase permuta- Figurc 52. The solutions to (217) and the distributions of initial guesses in the (132, Try-plane with three equal element lengths and 70:2 for material B. ............................................................................................................. ' ...... 1 17 Figure 53. The solution to (217) and the distributions of initial guessesin the xiii (32, Bag-plane with three equal element lengths and 70:08 for material B. ................................................................................................................... 118 Figure 54. The solution to (217) and the distributions of initial guesses in the (.112, E3)-plane with three equal element lengths and 70:08 for material A. All initial guesses flow downhill to the <1, I, I> solution under the it- eration scheme of Section 10.5. .................................................................... 119 Figure 55. The solutions to (217) and the possible initial guesses in the (172, Hap-plane with three equal element lengths and 70:2 for material A. ........................... 120 Figure 56. The solutions to (217) and the distributions of initial guesses in the (52, Kayplane with three equal element lengths and 750.97 for material B. The contour diagram for this case is given in Figure 50. ......................... 123 Figure 57. The permutationally distinct solutions for N from N=l to N=8 with ma- terial A and material B under the hard device loading. The number in- side the [.] indicates the number of permutationally distinct solutions. Note that the positions of Yoz are not in scale. The T give values of 70 at which the number of permutationally distinct solutions suddenly chang- es. Often the number of permutationally distinct solutions exactly at the T values is different from the number of permutationally distinct solu- tions that is available on either side of the T values. .................................... 128 Figure 58. The solution paths in 912 obeying the minimum energy solution criterion for the two-element, the four-element, the eight-element and the twenty- element cases under a monotonically increasing hard device loading with material A. ............................................................................................ 130 Figure 59. (a) The loading 0(t) given by (231). (b) The corresponding response function 70(t) which is given by (232) using (231) for material A with k=4.05. Three different total number of time steps, 10, 20 and 30, are ‘ considered in the numerical computation of this “exact” solution. .............. 137 xiv Fi Figi [:1ng FigUli Figure 60. The solution paths obeying the kinetic criterion for different algorithms under a hard device loading with 10 loading time steps. Here the condi- tions are given by (229) and (231) with material A and k=4.05. .................. 139 Figure 61. The effect of refining the spatial discretization. The solution paths obey- in g the kinetic criterion for the two-element (a), the four-element (b) and the eight-element cases(c). This is done by a hard device loading with three different loading time steps. Here the conditions are given by (229) and (231) with material A and k=4.05. The curves indicated by i are available for the value of N being considered. The combination algo- rithm was used for implementing the kinetic criterion. ................................ 141 Figure 62. The effect of refining the time discretization. The solution paths obeying kinetic criterion for the two-element, the four—element and the eight-el- ement cases. This is done by a hard device loading with three different loading time steps, 10 time steps (a), 20 time steps (b) and 30 time steps (c) for initial conditions given by (229) with material A and k=4.05. The combination algorithm was used for implementing the kinetic criteri- on. ................................................................................................................. 142 Figure 63. Three loadings with different loading rates. .................................................. 143 Figure 64. The effect of refuting the spatial discretization. The solution paths obey- in g the kinetic criterion for the two-element (a), the four-element (b) and the eight—element cases (c) under the three hard device ramp loadings. Each loading has a different loading rate for conditions given by (235) with material A and k=10. The curves indicated by l are available for the different element discretizations. The combination algorithm was used for implementing the kinetic criterion. ................................................. 144 Figure 65. The solution paths previously displayed in Figure 64 are here displayed (for y>l and material A). Here we group the FEM paths associated with XV the same loading rate on the same frame. Within each frame we see the effect of increasing the number of elements. ................................................ 145 Figure 66. The solution paths in (212 obeying the minimum energy solution criterion for the two-element, the four-element, the eight-element and the twenty- element cases under a monotonically increasing hard device loading with material B. ............................................................................................. 147 Figure 67. (a) The loading o(t) given by (231) which is also considered for material A. (b) The corresponding response function 70(t) which is given by (236) using (231) for material B with k=4.05. Three different total number of time steps, 10, 20 and 30, are considered in the numerical computation of this “exact” solution. ................................................................................ 149 Figure 68. The effect of refining the spatial discretization. The solution paths obey- in g the kinetic criterion for the two-element (a), the four-element (b) and the eight-element cases(c). This is done by a hard device loading with three different loadin g time steps. Here the conditions are given by (229) and (231) with material B and k=4.05. The curves indicated by l are available for the value of N being considered. The combination algo- rithm was used for implementing the kinetic criterion. ................................ 150 Figure 69. The effect of refining the time discretization. The solution paths obeying the kinetic criterion for the two—element, the four-element and the eight- element cases. This is done by a hard device loading with three different loading time steps, 10 time steps (a), 20 time steps (b) and 30 time steps (c) for initial conditions given by (229) with material B and k=4.05. The combination algorithm was used for implementing the kinetic criteri- on. ................................................................................................................. 151 Figure 70. The effect of refining the spatial discretization. The solution paths obey; ing the kinetic criterion for the two-element (a), the four-element (b) and xvi the eight-element cases (c) under the three hard device ramp loadings. Each loading has a different loading rate for conditions given by (235) with material B and k=10. The curves indicated by l are available for the different element discretization. The combination algorithm was used for implementing the kinetic criterion. ................................................. 153 Figure 71. The solution paths previously displayed in Figure 70 are here displayed (for y>1 and material B). Here we group the FEM paths associated with the same loading rate on the same frame. Within each frame we see the effect of increasing the number of elements. ................................................ 154 xvii 1. Introduction 1.1 Statement of problem Elasto-(quasi)-static fields associated with discontinuous strain in the fully nonlinear equi- librium theory of elasticity have been explored recently [1987A] in order to model phase transitions. Stress-induced phase transitions in solids can be modeled in a continuum elas- ticity setting by using non-convex strain energy density functions. For an elastic solid, if its strain energy density function is non-convex, then its displacement equation of equilibrium loses ellipticity and many new features of the theory emerge. For example, discontinuous strains can occur, solutions are less smooth than usual, and massive failure of uniqueness can occur in even the simplest boundary value problem. In addition, a breakdown of the familiar balance between the rate of work and the rate of increase of stored energy can oc- cur [197 9K]. In particular, the departure from the conventional elastic behavior arises in quasi-static motions of an elastic body in which, at each instant, the associated equilibrium field involves discontinuous strain. The purpose of this study is to investigate the connection between analytical re- sults and numerical results of phase transition behavior within a strain-softening bar sub- jected to either a given displacement or to a given load at ends. In particular, the bar may have a nonunique strain which is due to the strain-softening of the material. Moreover, we study loading behavior under two separate criteria for resolving the nonunique issue mentioned above, one is the minimum energy solution criterion, in which the bar always assumes configurations that have the minimum energy, the other is the kinetic criterion, in which the velocity of phase boundaries in the bar are determined by a kinetic relation. Finally we compare the analytical results and the numerical results that are obtained by using the finite element method. l.2 Literature review A material is said to be in a strain-softening state when the product of its strain rate and its stress rate is less than zero. Therefore the strain-softening takes place when the stiffness matrix is not positive definite. Strain-softening can be found in certain materials such as concrete, rocks, dense soils and some composites [1984R], [19850]. Shan et a1 studied the fracture behavior of a “quasi-brittle” cement-based material by experimental methods [19938]. The “quasi-brittle” material, which is different from a brittle material, responds nonlinearly prior to the peak stress and gradually decreases its stress after reaching the peak. One finds strain-softening in the constitutive equations for phenomena of damage models where the strain-softening is due to plastic fracturing [1979B]. Here the final col- lapse of compression test of concrete often happens after the maximum stress is reached. Both the size effect and the shape effect play important roles in the strain-softening of con— crete. Belytschko et al studied transient solution for one-dimensional problems with two types of stress-strain relations: one where the stress remains constant after the strain-soft- ening and the other where the stress increases again [1987B]. Bazant studied wave propa- gation within a strain-softening bar which is loaded by forcing both ends to move simulta- neously outward, with constant opposite velocities. The analytical solutions and numerical solutions which are obtained by using the finite element method are given [19853]. A simple model is useful since strain-softening involves issues related to material instability, nonuniqueness and existence. A one-dimensional bar is the simplest model studied for an elastic solid. Ericksen [1975E] considered a bar whose stress-strain rela- tion is non-monotonic and explored it for material instability by using the equilibrium theory. This can be viewed as a simple model for phase transition in an elastic solid. Non- monotonicity of the stress response behavior in [197 5B] may lead to three allowed values of strain which correspond to a single value of stress. These strains are in three different strain phases: low strain phase, intermediate strain phase and high strain phase. A prop- agating strain discontinuity in these material may be either a shock wave or a phase boundary, according to whether the locations separated by the discontinuity are of the same phase or different phases. In [1975E] the intermediate strain phase is a phase in a behavior of strain-softening. To formulate the problem, the formal presentation of basic laws is useful. Gurtin [19916] proposed formal presentation of fundamental laws in a fully dynamical setting: balance laws for mass, forces and energy and a law of the entropy growth for a two-phase continuum with sharp phase-interface. In addition, he explained the ideas involved in the formulation of these basic laws and clarified the differences between various formula- tions. The problem can be studied in either an equilibrium setting, a quasi-static setting or a dynamical setting. The quasi-static setting gives motion by advancing through a stat- ic case at each instant so that the difference between a quasi-static setting and a dynam- ical setting is the absence or presence of inertia. The dynamical setting usually involves wave propagation but there are no waves in the quasi-static setting. One might expect that the quasi-static setting would give the large time solution after waves have damped out. If body forces are absent and the field is in equilibrium, then the equilibrium equation can be expressed in terms of stress a as V0 0 = 0, where V0 is divergence op- erator, and o is the Piola stress which is the force per undeformed area [1977K]. Hence for an one-dimensional system the equilibrium equation can be expressed in terms of stress 0 as o’(X) = 0, here ’ denotes differentiation. The equilibrium equation may be locally elliptic, parabolic or hyperbolic depending on the material and the deformation. Knowles and Stemberg (see [1978K]) showed that the equation type was established from the local slope of the stress-strain relation as follows 6’ (y) > O 4:» elliptic, 6’ (y) = O aparabolic, 6’ (y) < O c: hyperbolic, where 6 (y) is stress and ‘y is strain. If the equilibrium field contains a curve along which the partial derivative of displacement fails to be continuous, the curve is known as an equi- librium shock or a stationary phase boundary. Equilibrium shocks affect the balance of the energy in the elastostatic field [1977K], [1978K] and [1979K]. Abeyaratne [1983A] has derived “equilibrium shocks” on which the deformation gradient and stress are discontinu- ous across certain internal material surfaces. Equilibrium equations, natural boundary con- ditions, traction continuity condition and a supplementary jump condition are also derived. Finally he concluded that this supplementary jump condition implied that a stable equilib- rium shock must necessarily be dissipation-free. Kikuchi and Triantafyllidis [1982K] also explored the admissible equilibrium solutions to the Euler-Lagrange equation for a mini- mization problem of the energy. Abeyaratne and Knowles [1987AA] have investigated a twist problem in the the- ory of finite deformations of elastic materials whose associated equation of equilibrium does not remain elliptic at all strains. Abeyaratne and Knowles [1989A] have concluded some elastic materials are capable of sustaining finite equilibrium deformations with dis- continuous strains. Boundary value problems for such “unstable” elastic materials often pose an infinite number of solutions. Abeyaratne and Knowles [1987A] have considered maximally dissipative quasi-static motions and found the possibility of modeling dissi- pative mechanical response in solids, on the basis of the equilibrium theory of finite elas- ticity for materials that may lose ellipticity at large strains. The result permits the occur- rence of discontinuous displacements gradients. For quasi-static motions, one can obtain an unique solution by requiring that all equilibrium states are minimizers of the potential energy functional. The present value of twist is then uniquely determined by the‘present value of torque and vice-versa, as in conventional elasticity. But in contrast, for maxi- mally dissipative motions, twist depends on the past history of torque during the motion, as in plasticity. Abeyaratne and Knowles [1987AA] have concluded that quasi-static motions, for which the equation of equilibrium does not remain elliptic at all strains, may be dissipa- tive. In particular, they show that the macroscopic response in the problem treated may correspond to elastic-perfectly plastic behavior. It was shown that the quasi-static re- sponse of a bar to a prescribed force history is then fully determined if one considers ad- ditional information in the form of a kinetic relation [1989A]. In particular, they show how unstable elastic materials can be used to model macroscopic behavior as in visco- plasticity. Neither the equilibrium setting nor the quasi-static setting treats truly dynamical problems involving wave propagation. Due to nonuniqueness of dynamical solutions, a number of authors have studied many criteria to single out an unique solution. All entropy is defined to be a ratio between the heat transfer to the system and. the absolute tempera- ture in a reversible process. Moreover a dissipation, or degradation, is defined to be the difference between the maximum possible work and the actual output. In an adiabatic system the dissipation is proportional to the increase of entropy. In this setting the prin- ciple of increase of the entropy is simply another statement of dissipation of energy. Both of which are also equivalent to the second law of thermodynamics. Dafermos [1973D] extended the entropy admissibility criterion of Lax [1973L] to entropy rate admissibility criterion for solutions of a hyperbolic system of conservation laws. The selected solu- tions are those in which the total entropy decreases with the highest possible rate. Dynamical phase transitions are also a highly studied topic in gas dynamics [19831-1], [19861-1], [1991AAA]. The van der Waals fluid problem could be viewed as a counterpart of the phase transformation problem in finite elasticity. Moving phase boundaries are kinematically similar to other types of shock fronts in that velocity and strain are discontinuous across them. Hagan and Slemrod [1983H] have investigated ad- missibility criteria for nonlinear conservation laws based on capillarity and viscosity, and found results consistent with experiment for materials which exhibit phase transitions by means of specific examples. Here an important concept is the Maxwell stress 6 = o’maxwell’ which is uniquely defined by the requirement that the two loops enclosed by the stress-strain curve a = 6 (y) and the horizontal line 0 = 0' we" are equal in area. Hattori [19861-1] has treated special Riemann problems and shown that the well known Maxwell construction is admissible, according to the entropy rate criterion. Af- fouf and Caflisch [1991AAA] have presented a numerical study of isothermal fluid equa- tions with a nonmonotone equation of state. It can serve as a model for describing dy- namic phase transitions, and the model can explore an analytic argument for stability of phase transitions. The analysis of the fully dynamical motion of a phase boundary for a specific class of elastic materials whose stress-strain relation in simple shear is nonmonotone can be found in Pence [1991P], [1991PP]. It was shown that the ringing of a shear pulse be- tween the external boundaries and an internal phase boundary gives rise to periodic phase boundary motion for both the case of a completely reflecting phase boundary and the case of a completely transmitting phase boundary [1991P]. The possibly infinite solutions that are not completely reflecting or completely transmitting are examined from the perspec- tive of energy and dissipation. It is shown that there also exist at most two solutions which involve no dissipation, and one solution that maximizes the mechanical energy dissipation rate [1991PP]. Similar topics in one-dimensional elasticity can be found in [1980]], [1986P], [1987?]. Pence investigated emergence and propagation of a phase boundary in an elastic bar, whose stress-strain relation is nonmonotonic [1986P]. An asymptotic description of the emergence of the phase boundary was provided for short times. The energetic behav- ior of a phase boundary that is subjected to concurrent dynamic pulses under a maximally dissipative criterion and under a kinetic criterion is studied by Lin and Pence ([1992L], [1993L]).The results show that the total energy loss is different from that which would occur if the two pulses were not concurrent. The dissipation which is due to wave ringing in non-elliptic elastic materials for large time is explored by Lin and Pence ([1993LL]. They set up a framework for treating the resulting wave reverberations and computing the energy dissipation for large time. Their numerical results shew that the large time en- ergy dissipation matches that which is necessary to settle into a new energy minimal equilibrium state. They obtain the same results analytically for the special case of a small dynamical perturbation. A very elegant mathematical investigation of a related problem can be found in Ball et al [1991B]. They explore mathematical models for the dynamical behavior of me- chanical systems that dissipate energy as time increases, and establish global existence and uniqueness results for the long time behavior of the systems. This provides insight into the dynamical development of finer and finer microstructure by certain material phase transformations. The lack of uniqueness of solutions to these problems can also be resolved by im- posing two constitutive requirements, a nucleation criterion and a kinetic relation [1988AA], [1990A], [1991A], [1991AA]. A kinetic relation can be in a form f = (p(s'), where f is the driving traction acting at the phase boundary and s’ is the phase boundary velocity. The relation between the driving traction on the phase boundary and the phase boundary velocity is similar to the relation between maximum dissipation to the maximum plastic work [1992A]. If one adds the viscosity and second strain gradient to the elastic part of the stress, it is equivalent to the imposition of a particular kinetic relation at the phase boundaries that propagate subsonically [1988AA], [1991A]. In a general thermo- mechanical process, the kinetic relation could be function of the driving traction and tem- io d1 .Bal chal and into Phat P031 1198 “her VCIOQ boun lflasu panC Dhas6 niech ior of a phase boundary that is subjected to concurrent dynamic pulses under a maximally dissipative criterion and under a kinetic criterion is studied by Lin and Pence ([1992L], [1993L]).The results show that the total energy loss is different from that which would occur if the two pulses were not concurrent. The dissipation which is due to wave ringing in non-elliptic elastic materials for large time is explored by Lin and Pence ([1993LL]. They set up a framework for treating the resulting wave reverberations and computing the energy dissipation for large time. Their numerical results shew that the large time en- ergy dissipation matches that which is necessary to settle into a new energy minimal equilibrium state. They obtain the same results analytically for the special case of a small dynamical perturbation. A very elegant mathematical investigation of a related problem can be found in Ball et al [1991B]. They explore mathematical models for the dynamical behavior of me- chanical systems that dissipate energy as time increases, and establish global existence and uniqueness results for the long time behavior of the systems. This provides insight into the dynamical development of finer and finer microstructure by certain material phase transformations. The lack of uniqueness of solutions to these problems can also be resolved by im- posing two constitutive requirements, a nucleation criterion and a kinetic relation [1988AA], [1990A], [1991A], [1991AA]. A kinetic relation can be in a form f = tp(s'), where f is the driving traction acting at the phase boundary and s is the phase boundary velocity. The relation between the driving traction on the phase boundary and the phase boundary velocity is similar to the relation between maximum dissipation to the maximum plastic work [1992A]. If one adds the viscosity and second strain gradient to the elastic part of the stress, it is equivalent to the imposition of a particular kinetic relation at the phase boundaries that propagate subsonically [1988AA], [1991A]. In a general thermo- mechanical process, the kinetic relation could be function of the driving traction and tem- perature [1990A]. It is also possible that a modified version of the entropy rate admissi- bility criterion can act as a particular kind of kinetic relation. For complicated boundary value problems, it is very difficult to find analytical so- lutions. Therefore we need numerical analysis to obtain approximate solutions. Numeri- cal analysis using these types of phase transition models can be found in [198888]. Sill- in g applied numerical analysis to plane strain equilibrium problems in finite elasticity us- ing a finite-difference dynamic relaxation method. He has claimed this method is useful for the study of localization and phase changes in elastic solids. Silling [19888] has also applied numerical analysis to the equilibrium anti-plane shear deformations for materials in which the governing equation varies in type locally from elliptic to hyperbolic. This extends Ericksen’s analysis for bars to the two-dimen- sional case by numerical analysis. He has shown the emergence of surfaces of disconti- nuity in the displacement field in a weakening material and in a trilinear material, and has also found the solution to equilibrium problems by using the dynamic relaxation method. 1.3 Section overview It is obvious from the literature review that nonuniqueness of the strain-softening material is an important issue which is not yet fully resolved. In order to better understand the non- uniqueness of the strain-softening material, we choose a simple system, a bar, to model the behavior. We study the phase transition behavior of the bar both analytically and numeri- cally. Unfortunately, there are few numerical analysis results obtained by using finite ele- ment method (FEM) for this kind of problem. One of the purposes of this thesis is to study the connections between the analytical results and numerical results which are obtained by using FEM. We are concerned with the phase transformation induced by stress within a strain- softening bar. The main goal is comparing the solutions of the analytical study with those of the numerical study. That is why we devote nine sections, from Section 1 to Section 9, to the analytical study. This dissertation is divided into sixteen sections. The first section is an introduc- tion. The next fourteen sections, from Section 2 to Section 15, are grouped into two parts. In Part I, from Section 2 to Section 9, there is an analytical study of a continuous system. In Part II, from Section 10 to Section 15, there is a numerical study (FEM) of the associ- ated discrete system. Two example materials: material A and material B, meaning two specific nonlinear elastic response functions, are considered in detail. The material prop- erties which give strain-softening behavior are described in Section 2 (see Figure l). Phase-I is defined to be the phase of strain between 0 and 1, while phase-II is defined to be the phase with strains greater than 1. Section 2 derives the elastic stress response den- sity function, the corresponding strain energy function, and an elastic secant modulus function. In Section 3, two loading devices, the soft device and the hard device, are intro- duced. For the soft device the bar is firmly fixed at one end while the other end is sub- jected to a given dead load. Correspondingly, for the hard device the bar is also firmly fixed at one end while the other end is subjected to a given displacement. If we consider the material described in Section 2, then at each point in the bar there are two possible phases, either the lower strain phase (phase-I) or the high strain phase (phase-H), to choose from. For a bar having, at most, two regions in different phas- es, there exist four possible ways to arrange these phases: (1,1), (1,2), (2,1) and (2,2). For example the (l,2)-case is a state with the lower strain phase on the left and the high strain phase on the right. Section 4 is devoted to the study of two phase equilibrium configura- tions with a single phase boundary. The energy functions for the soft device and hard de- vice are also derived. 10 In Section 5 we introduce and discuss the solution region for problems with a sin- gle phase boundary. This is done for both devices. By solution region we mean, in brief, a parameter plane region of stress 6 vs. average strain 70 for possible solutions. We study the solution region in the (1,2)-case. In particular we first explore the lines of constant stress 0 where o is defined in (3). Secondly we explore the lines of constant strain 70 where 70 is defined in (18). Finally we explore the lines of constant s where s is a non- dimensional location parameter for the phase boundary and is defined in (22). The corresponding hard device energy functions on the boundary, lines of constant 0, lines of constant 7, and lines of constant 5 within the solution region in the (1,2)-case are discussed in Section 6. The minimum energy criterion is introduced at the end of the section. Similarly, Section 7 is devoted to the soft device. Quasi-static motion is defined and introduced in Section 8. The dissipation and the driving traction at the phase boundary of a bar are discussed. We study the admissible direction of movement of the phase boundary by the second law of thermodynamics un- der isothermal conditions. Finally we study the admissible paths which have a non-neg- ative product of the phase boundary velocity and the driving traction at the phase bound- ary. Recall from literature review that if the local slope of stress-strain relation is neg- ative then the equilibrium equation is non-elliptic. For a non-elliptic equilibrium equa- tion the solutions are not always unique. Since the massive failure of uniqueness is relat- ed to the stress-strain relation which is the constitutive relation of material, this kind of non-uniqueness has been viewed by some as a constitutive deficiency. To resolve a con- stitutive deficiency, one may supplement the theory with a kinetic criterion which applies to phase boundaries only [1989A]. The kinematics of the phase boundary and the. admis- sible paths with such a kinetic criterion are discussed in Section 9. 11 To analyze this boundary value problem (sketched in Figure 3) numerically we choose the finite element method. At first the basic finite element formulation of a bar element is derived in Section 10. Second, we select two-node elements to compute the interpolation functions (160) and the element equation (163). Third we assemble the el- ement equations of each element into a resulting system equation (175). After this we consider boundary conditions for both the hard and soft devices to find a reduced system equation ((184) or (188)) for a general N ~element model. Finally we discuss an iterative scheme to find the solution from an initial guess. In Section 11 we first find the solutions of the system equation for a two-element model. We study the relations between the solutions and the possible nodal displace- ments to explore the nonuniqueness. This simplest of all possible cases allows us to ob- serve how the nonuniqueness in analytical solutions translates into nonuniqueness in nu- merical FEM solutions. Then we investigate the corresponding energy as a function of the possible nodal displacements and relate this to the convergence and divergence of the iteration scheme. Similarly in Section 12 the available solutions for a three-element model are dis- cussed as in Section 11. In Section 13 the available solutions for N >3 are discussed. For resolving the nonuniqueness issues mentioned above we explore the numeri- cal results using the finite element method under either a minimum energy solution cri- terion or else a kinetic criterion. Several algorithms are chosen to implement the kinetic criterion. Then the results obtained by using these algorithms are compared. This is done for material A in Section 14 and for material B in Section 15. In particular, for the kinetic criterion we compare the exact solutions and the solutions which are obtained by using the finite element method. Finally we determine the FEM solutions obtained using the kinetic criterion for three ramp loadings with different loading rates. 12 In Section 16 we have the conclusions from this study and recommendations for future work. 2. The material model In this thesis we consider a compressible homogeneous isotropic hyperelastic bar with unit cross section area. The bar in the reference configuration is described by the region R = {0 s X S L} . In the deformed configuration the particle at X in the undeformed state is moved to X + u(X), where u(X) is the displacement of the particle originally at X. The strain ‘7 = du(X)/dX of the particle at X is restricted to the range ‘7 2 —1 in order to assure that the mapping X —> X + u(X) is invertible; in fact, we shall limit attention to 72 O. The elastic stress response function of the bar is taken to be in the following: ‘ = E“, s , DC?) = 1(7)“ 7 A? 7M (1) 611(7), Y>7Mt where 6,,(?) will be given in what follows. Note that 6, (y) is a linear function of ii and E is a positive constant, but 611(7) may not be. In particular, the behavior of a strain-soft- ening 6” (’y) , a behavior in which the stress declines for increasing strain, is investigated. Note fiom Figure 1 that "o I I (’j') is strictly decreasing (and hence invertible) and has the fol- e5M ‘0 (i) a] (it) 3,, (7) 9c .................................................................................................. ¢ > Figure l. The stress response function of the bar with positive strain. 13 14 lowing properties: be < 6” (’9) 5 6M, ’6”’ (y) < 0, 6”” (y) > 0, 6” (y) —>6c20, when ’y—Mo. (2) Here ' denotes differentiation. The stress response function of the bar is described by two strain phases: the phase of strain from O to if)" is called the lower strain phase (phase-I) and the phase of strain greater than ’iM is called the high strain phase (phase-II). If two non-dimensional parameters are introduced _ fr . 0:. . (3 M .g ll i" I <> Q then the non-dimensionalized stress response function (see Figure 2) of the bar may be re— writtenas o = _ 7 011(7), y> 1. A l 6 _ 5() (Y) 01”) ll 7 o ................................................................................................... C < J Figure 2. The non-dimensionalized stress response function of the bar which is used in this thesis. 15 Note that 61(7) is also a linear function of y and 611(7) has the following properties: oc< 611(7) 51, 61]’(Y) < 0! 6””(7) > 0, 611(7) —->oc20 when y—->oo, (5) where CC = Ere/33M. Now phase-I is the phase of strain from O to 1 and phase-II is the phase of strain greater than 1. Since both 61(7) and 611(7) are invertible, the associated inverse functions can be rewritten in terms of o. This gives the corresponding inverse functions 7 = rim), i = I or II. (6) Note from (4)1 that F,(o) = o. The corresponding normalized strain energy density func- tion W(y) is given by W(Y) = 1 (7) "2’+I‘1!6[](‘Y)dyr 7)], so that the stress is equal to the first derivative of the strain energy density function, 5(7) = W'(Y)- (8) Note since 6,,(7) —> at 2 0 as y—> co that lim WW) > 0. (9) t -) °° however the limit lim (W (7) — 706) a We > O (10) r -* ~ may or may not be finite depending on the rate at which 611(7) —> cc. Here We is the area under the stress-strain curve which is also above the horizontal line a = O'c. Two particular functions for 611(7) are given in Appendix-A and Appendix-B. 16 These functions are used from time to time in examples in this thesis. Thus we say that these functions define two example materials: material A and material B. Both of these materials have CC = 1/ 2, however material A has WC = oo while material B has a finite value of WC. Let 000 be the stress field in the bar, then the equilibrium equation, with no body force, is given by —d—o(X) = 0. (11) dX Therefore, the stresses along the bar are spatially constant, 6(X) = constant, 0 SK SL. (12) In what follows it will be convenient to define an elastic secant modulus function by 5(7) = 5(7)”. (13) giving E,(y) = 1, ygl, 5(7) = 5 (Y) (14) E,,(y) = —”—, y>1. 3. Two loading devices In this study two standard types of loading procedures, the soft device and the hard device as depicted pictorially in Figure 3, are considered. 3.1 Soft device (Figure 3) Here the bar is firmly fixed at one end while the other end is subjected to a given dead con- stant load F = F o = A00 2 0’0. Since the cross sectional area A = l, the boundary con- ditions are u(O) = 0 and (Six = L = 00. This is the so-called force problem that consists of finding an unknown displacement field u satisfying the boundary conditions. The solution is not unique in this problem. For example: 00 < 0c =9 1 solution, cc < 00 < 1 => co solutions, (15) 00 > 1 = no solution. / X ; u(O) = 0 F = 0,u(L) = 0 undeformed 7 ]--> 7 a Z “(0) = 0 F = F0,u(L) = ? deformed: 7 1., soft device 7 a = F = ?,u(L) = 80 deformed: g “(0) 0 hard device a /7 Figure 3. Undeforrned and deformed configurations of soft device and hard device. 17 18 From (4), (5) and Figure 2 the only admissible value of strain whenever 0’0 < do is y = 0'0, so that the bar has one solution which is phase-I everywhere. Hence 0C has the interpreta- tion of being the low limit on the stress that can support phase-II. It is impossible to have solutions for o> 1 because the bar cannot resist the stress greater than 1. The case ac < o < 1 is more complicated and interesting. The bar could be in a mixed state of phase- 1 and phase-II, since then two values of strain are compatible with the one value of stress. The final elongation will then depend on the fraction of the bar which is in phase-I versus the fraction of the bar which is in phase-II. In particular, if the bar is partitioned into n sec- tions, then each section could be in either phase-l or phase-II. In fact, there would be 2" choices. Therefore, if the number of sections is increased, then the number of choices is in- creased too. In this way an infinite number of solutions could be constructed. For the force problem, classical absolutely stable equilibrium configurations must minimize the follow- ing soft device energy function can (L) L a = wow-’5‘) - (16> OT—‘h‘ 3.2 Hard device (Figure 3) On the other hand, if the bar is firmly fixed at one end while the other end is subjected to a given constant displacement 80 > O , then the boundary conditions are u(O) = 0 and u(L) = 50. Let 70 E So/L be the average strain in the bar. This is the so-called elongation problem that consists of finding of an unknown constant force F = 60 and a displacement field u(X) satisfying the boundary conditions. Again the solution is not unique in this prob- lem and it will be described in more detail in Section 5. For this problem, classical abso- lutely stable equilibrium configurations must minimize the following hard device energy function l9 1 A = ([Wemif). (17) In what follows the applied load 60 will be replaced by 0' for simplicity. In this study, the hard device problem is our main focus, and unless it is to remark otherwise, the reader should assume that the hard device case is under consideration. 4. Two phase equilibrium configurations with a single phase boundary The nonmonotone stress response function permits two values of strain to give the same value of stress. Thus solutions (12) to (1 1) may involve alternating regions along the bar of these two strain values. lrnagine that there are two different strain regions in the bar, each with constant strain 7L and 7R separated at X = h with 0 < h i = I, j = II, 21 (2,1)-case when 7L2] and OSYRS 1 =>i = II, j = I, ( ) (2,2)-case when 7L 2 1 and 7R 2 1 =>i = 11, j = II. For the (l,l)-case and the (22)-case, the strains 7L and 7,, are Qua] and so the location h is not of significance. However it will be important to consider these two “trivial” cases in what follows as naturally limiting states when either h -> 0 or h —> L. The location X = h in the (1,2)-case and the (2,1)-case have the interpretation of being a phase boundary. It is useful to introduce a non-dimensional phase boundary location parameter s = E (22) then by virtue of the boundary conditions u(O) = 0 and u(L) = 80, it can be obtained that 20 21 8 rLs+rR(1-s) = 1357,. (23) Now, by using (17), the hard device energy function of the bar is given by s 1 X X A, = jW,(YL)d(Z) + jwmmzl y: ray orzan, 0 s 24 i=1(1)or2(11). ( ) = Wms + Wm) (1 - s) = 3,,(0. 7,) In this thesis, a superposed A will indicate that the quantity is regarded as a function of o and 70. For AU“), 70) this requires the use of (6) for 7L and 7R; also the elimination of 5 uses (23) to define the function 3,1(6, 70), where _ YR - 70 13(0) " Yo s " 'YR - Yr = 13(0) - Fla) E 3"” 70) ' (25) Similarly, from (16), the soft device energy function of the bar is provided by .t 1 5,, = {la-(mag) + [ijpdt’l—f) — or, ,-= 1(1) 0, 2(11), (26) S A j = 1(1) or 201). = wings + W,(YR) (1 - s) - or, = ~=r,(°t 7.) Note that E. {1(6, 70) is also a function of o and ya. 5. The solution region with a single phase boundary The two-phase hard device problem involves finding a satisfying (23), (6) and the condi- tion 6(7L) = 6(7),). Here the freedom to let s take on any value in 0 < s < 1 gives rise to similar nonuniqueness phenomena as discussed for the soft device problem in Section 3.1. In particular, given 70, there is often a range of a that satisfy the above conditions. This range shall be denoted by [1(70). As yo increases from 0, it can be imagined that this set l'l(yo) sweeps out a solution region Q = [ (yo, 11(7)) l 70 > 0} in a (70, 0') -plane. This “solution region” representation will be extremely useful in this thesis and so will be de- scribed now in some detail for each of the four cases given in (21). The solution regions will be denoted by Oil. for each separate case. Solutions only exist for the (l,1)-case if 0 < yo < 1. Then the range 11(70) collapses to a point and the associated solution region is the straight line (see Figure 4) £211: {(yo,o)|o=6,(yo),0 1. Then again the range 11(70) collapses to a point and the associated solution region is the curve (see Figure 5) 022 = {(70, 0)] o = 6,,(70), yo >1}. (28) On the other hand for the (l ,2)—case, solutions exist if 70 > 0c and are given by (see Figure 6): Figure 4. The solution region Q”. 22 23 Figure 5. The solution region 922. £212 = {(70, 0)] oc< 0< 6,(yo),if 0c 1. Foragiven (yo, 0) e (212, the associated phase boundary location s is given by (25), and values of s on the boundary of 912 are shown in Figure 6. Similarly, it is found that (2:9 21 (30) 12‘ However these solutions involve 521 —) l - s12 so that the values of s on the boundary of (221 are shown in Flgure 7. Finally the solution regions for the all four cases are shown in Figme 6. The solution region 912. 24 e A 1 ore/é KW") a, WIH;.///////%.%//S/ 1 . ’g 1: 70 > Figure 8. The solution regions (2", 9,2, (In and (In. Figure 8. The above discussions of Q Q 021 and 922 are based on the hard device 11’ 12’ problem, however, they also can apply to the soft device problem. That is, given 0, there is also a range of 70 that satisfy (23). And according to 0 < s < 1 the value of 70 is always between I‘,(0) and I‘”(0). As 0 increases fromO to 0c there is an unique solution 0 = 10. While as 0 increases from 0 c, it is easy to imagine the range of 70 which satisfies (23) sweeping out the same solution region. In what follows, the (12)-case is studied. The following discussion will further 25 make clear why (212 has the form given in Figure 6. 5.1 The boundary of £212 From (29) and Figure 6 the boundary of £212 consists of three curves: 0 = 01(70), 0 = 6,,(70) and 0 = 0C. These three separate curves can be considered as limiting cases and are described in the following: The curve 0 = 0,010) coincides with Q“ . Therefore, it corresponds to a case in which the phase-II region has vanished. Since this phase-H region is to the right of the phase boundary for the (l,2)-case, the phase boundary has squeezed the phase-II region all the way to the end X = L, thus s = l. The solutions along the curve = 6,0) are therefore pure phase-1. Let’s define $1 to be this set of ordered pairs (70, 0) which are pure phase-I so- lutions on the interval 06 < 0 S 1: Q1: {(70.0)| (10:0,0C<051)}. (31) Similarly, the curve 0 = 61,00) coincides with (222. So it corresponds to a case in which the phase-I region has disappeared. Since this phase-I region is to the left of the phase boundary for the (12)-case, the phase boundary has squeezed phase-I region all the way to the end X = 0, thus s = 0. The solutions along the curve 0 = 0,,(70) are pure phase-II solutions. On the other hand, the curve 0 = 0 c has the interpretation of being an unattainable limit corresponding to an infinitely strained phase-II region of vanishing length. Like (211 it also corresponds to a case in which the phase-II region has vanished, but it doesn’t coin- cide with the pure phase-I solutions of (In. To understand this, recall that the total elonga- tion is given by (23). Since phase-II is to the right of the phase boundary, the phase. bound- ary squeezes the phase-II region all the way to the end X = L when s —-> 1. But the main 26 difference from Q“ is that the phase-I region on the left (0 SX _<. h) has the elongation 5:) = F,(0 9311th = 0 CL while the phase-II region on the right (h s X S L) has the elon- . r _ . . - . = 00 o _ = gauon 50 - olgnocf‘lfl0)hll_l)nL(L h). Here chinocl‘um) and Itth (L h) 0 so that the product which defines 5; corresponds to 00 ~ 0 where the rate of the one limit with respect to the other is chosen so that the product 5; remains fixed. Moreover the total elon- . . I r . . gation at L 15 matched to 50 + 80 = 50 = 70L. For different pornts (yo, 0) along the curve 0 = 06 the elongations 5: are always fixed at 5:) = 06L while the elongations 8; take on different values to match their corresponding 70, that is 6; = VOL—5; = (70 - 06) L which is greater than zero since these solutions only exist for 70 > 06. 5.2 Lines of constant 70 in (2,2 The lines of constant 70 in 912 imply that the elongation is prescribed while there is a range of nonunique 0 which satisfy (25). Prescribing 70 is the definition of the hard device. Therefore, the lines of constant ya in $212 represent solutions to the hard device problem for given average strain 70 = 5 0/ L associated with elongation 50. If the bar is elongated with 50 < ‘7ch then the only possible solution corresponds to on . This means the elonga- tion is not enough to reach phase-II. As a result, the bar is in pure phase-I and there is no phase transition at this stage. If the bar is elongated with 0CL < 80 < L, then the bar could still be in pure phase-I. However it could also be in a mixed state of phase-I and phase-II, but it cannot be in pure phase-H. On the other hand, if the bar is elongated with 80 > L, then the bar could be in a pure phase-11 solution; it could also be in a mixed state of phase-I and phase-II but it cannot be in pure phase-I. 27 5.3 Lines of constant 0 in (2,2 The lines of constant stress 0 with 0 > 0c in $212 represent possible solutions to the soft device problem for given F = 0 > 06. If the force F = 0 is prescribed with 0 < 06, then the only possible solution corresponds to 911. This also means that the prescribed force is not strong enough to reach the phase-II. As a result, the bar is in pure phase-I and there is also no phase transition at this stage. But if the bar is prescribed by the force F = 0 with 06 < 0 < 1, then there are an infinite number of possibilities corresponding to points in the lines of constant 0 in (212. The bar could then be in pure phase-I, pure phase-II or a mixed state of phase-I and phase-II. 5.4 Curves of constant s in (2,2 It is to be noted from (25) that _ .. _ r11(°)‘70 S — 312(69 YO) " rll(o)_rl(o.) (32) is a function of ya and 0 in the (yo, 0) -plane and it is easy to imagine the level curves of fixed s in $212. In this section the form of these curves in (212 is investigated. An example of the level curves is given in Figure 9 for a particular material (this particular choice of the function 6,,(7) is given by (a. l) of Appendix-A). It will be convenient for our subsequent discussion to reorient to solution region (212 so that stress 0 is the horizontal axis and av- erage strain 70 is the vertical axis. The reorientation that this gives for the constant s curves of Figure 9 are shown in Figure 10. For each fixed value of s e (0, 1) , let ya = 70(0) give the solutions to (32) so that the function ?0(0) defines a smooth curve on the (0, yo) -plane. There is a family of such curves on the (0, 70) -plane. Each curve corresponds to a differ- ent value of s e (0, 1) and it can be shown to be a single semi-infinite curve. The endpoint 28 O = 611(70) s=0 0.75 " f7 ._ s=l 0 = 61(7) r s: 19/20 0.5:0’ _ ........................................ Figure 9. The curves of constant s=(0.05, 0.1, 0.15,..., 0. 95) in the solution region (212 with the choice of the function 0,,(‘y_) given by (a. 1) of Appendix-A. The lim- iting case of s=0 corresponds to 0:0,,(70) and the limiting case of s=1 corre— sponds to both 0='y, with “>0: and 0:0,. 29 l r l 0 0.25 0.5: O": 0.75 O 1 Figure 10. It is convenient to reorient the solution region of Figure 4—Figure 9 so that the constant s curves display a minima. We shall use this orientation for the solu- tion region (212 in the rest of this thesis. 30 of each curve is the point at (0, yo) = (1, 1) , this curve is unbounded for large 7 and in the limit as y —> oo it tends to the asymptote 0 = 0 C. Equation (23) gives that 370 .5; O = yL—yR, on se (0,1). (33) Since 7L 5 1 _<. 7,, it follows that the different members of this family of curves do not in- tersect each other except possibly if 7L = 7R => 7L = 7R = 1. Furthermore, a curve corre- sponding to a larger value of s lies under that of the smaller value of s. Although the ex- ample in Figure 9 is for a particular 011(7) with 0 c = l/ 2, similar properties of these con- stant s curves are found for other functions 6,,(7) of the type obeying (5). It is to be noted from Figure 10 that if s is sufficiently small, then the constant s curves are monotonically decreasing from (0, yo) = (0c,oo) to (0,70) = (1,1). However if s is sufficiently large, then these curves first decrease and then increase. There is a transition value of s, say s separating these two groups of s curves (for the exam- ITOM’ ple material in Figure 10 it is found that s = 2/ 3 ). For s > s the constant s curves trans trans ’ have an internal minimum point while for s S s they don’t. It will be useful to know trans the location of these minimum points to the curves of constant 5. Taking partial derivative with respect to 0 with fixed s in (23) gives 37 330 -.- ?o'm) = I‘,’(0)s+l‘”’(0)(1 —s), (34) and setting this equal to zero gives l‘,’(0)s+1‘”’(0)(1 -s) = o, (35) Of 1“111") = FUTO) _ r], (0) E s (0), ~ (36) S 31 BY as a condition for 5?; = 0. It is to be noted that the transition value of s is given by 3 SW. = s"(1). (37) To verify that it is minimum at (0, yo) satisfying (35) we take partial derivative with re- spect to 0 in (34) again. As a result, it is found that 2 5:770 .—. ?0"(0) = 1‘,"(0)s+l‘”"(0) (l —s). (38) Since 01(1) is a linear function and F,(0) is the inverse function of 0,0), I",(0) is a linear function too and 6117) = 1’ r]’(o) = 19 6]”(7) = 1",”(0) = O (39) Thus the combination of (39) and (38) gives = 1"””(0)(1 —s). (40) S 2 30.2 Yo 2 Since 1 — s > 0 it follows that aa—zyo and 1“,,”(0) have the same sign provided that 3 ¢ 1. 0 Note that I‘"(0) is the inverse function of 0,,(7) so that o = 6,,(7) = 6,,(1‘,,(o)). (41) Taking first derivative with respect to 0 from (41) gives 1 = 0,,’(y)1‘”’(0). (42) The combination of (5)2 and (42) now gives 611’”) < O, I‘ll’(0) < O. (43) Taking second derivative with respect to 0 from (41) gives 0 = 6,,”(7) (I‘,,’(o)) 2 + 6,,’(r)1‘,,”(o). (44) Since, by (5), 0,,’(y) < 0 and 6, {(7) > 0, it is found that 32 FI,”(0) > 0 (45) and so (40) yields 535—:70 > O. (46) 3 Thus (40) shows that the curvature restriction 6””(7) > 0 given in (5) ensure that there are no maximum points to the constant s curves. These curves monotonically decay to the min- imum point from (0, yo) = ( 0 c, co) after which they monotonically increase to the point (0, yo) = (1, l) . Note that if the curvature restriction 6””(7) > O in (5) is dropped then other interesting possibilities of solutions for the constant s curves will occur, but these will not be discussed here. The minimum points along the curves of constant s could be obtained by (23) and (36), namely _ l‘,(0)l‘”’(0) — 1‘”(0)l‘,’(0) 70': - I‘ll’(0) - 17(0) (47) Let’s define $2 to be set of ordered pairs (0, 70) which are minimum points for the vari- ous curves of constant s (see Figure 1 1). Then $2 is also a curve on the interval 0c < 0 S 1. From (47) it is found that I‘ 0)I‘ ’(0)—F 0)l" ’(0) g2 = “(5’70)”; K 1"’1’1'(0)-I‘III'((0) I ’OC0 and if 0c< 0< 1 then F,,(0) -I‘,(0)>0 and I‘,’(0) = 1 , this gives 33 70 = l‘,(0) = 0 0 0.25 0.5: O": 0.75 O 1 Figure 11. The $2 curve in the solution region 9.; (again the choice of the function 3,,(7) is given by (al) of Appendix-A). 34 dY —0 >0. (50) do a 2 This shows that 70 is increasing along $2 with increasing stress. It is to be noted that the point (0, 70) = (1, 1) is the right end point of G2 and thus gives the largest stress 0 on 92. Therefore, it is found that 70 S 1 along the curve $2. Since $2 is monotonically in- creasing we can invert (50) and so can parametrize this curve by 70, that is 0 = 092(70). In addition (36) is invertible and the inverse function shall be denoted 0 = 0¢2(s). 6. The hard device energy function on (2,2 The distribution of the hard device energy function in $212 can be described by (24), A12 = 1112(0. yo) = w,(r,(o))i,2(o, 70) + W”(1‘”(0)) (1 - sum, 70)) . Level curves of fixed A12 in the region 012 are functions of ya and 0. An example of these level curves is given in Figure 12 for a particular material, that is a particular choice of the function 0,,(7) (in fact it is the same material used for the constant s curves in Figure 9). For this material and each fixed value of A12 6 (0.125, co) , the equation A12(0, 70) = A12 defines a smooth curve on the (0, yo) -plane again. There is a family of such curves on the (0, yo) - plane. Each curve corresponds to a different value of A12 6 (0.125, co) . In what follows the hard device energy functions A12(0, ya) on the boundary of £212, on the lines of con- stant 0, on the lines of constant 70 and on the lines of constant s are examined. Note that the general discussion applies to materials for which the function 0,,(7) has the properties as given by (5). 6.1 The hard device energy function on the boundary of 912 First of all, the boundary yo = F,(0) of the solution region (212 is in phase-I. As a result, the hard device energy function will be A12 = W,(F,(0)). Secondly, the boundary yo = 1"”(0) of the solution region (212 is in phase-II. As a consequence, the hard device energy function will be A12 = W,,(I‘”(0)). Finally, recall that the boundary 0 == 0 c of the solution region {212 is the non-attainable limit corresponding to an infinitely strained phase-II region of vanishing length. It is to be noted that 71 = 0 C and 72 —> oo on the TR - 70 7R '71. =1. boundary 0 = 0 C in the (1,2)-case, and also from (25) it is verified that s = 35 36 3 l . I \An = 2.0 A = 13 12 70 . A12 — 19 A12 = 12 A12 = 18 A — 17 12 ° A12 = 1.1 A12 — 1.6 A12 - 10 A12 _ 15 — A — 14 2 _ A12 — 09 12 _ A12 — 0.8 A12 — 0.7 A12 = 0.6 A12 = 0.5 A12 = 0.4 1 - ., _ _ A12 = 0.3 A12 = 0.2 A12 = 0.125 o l l l o 0.25' 0.5: a 0.75 0’ 1 C Figure 12. The level curves of constant hard device egergy function A12 in the solution region (212 with the choice of the function 0,,( 7) given by (a. 1) of Appendix- A. 37 We now show that the hard device energy function A12 on 0 = 0c is given by A12 = 06 (yo — 0C/ 2) : from (24), (25) and (7) the hard device energy function is given by A12=‘5(Wl(rl( »,i—Y°,:+Wu 7A . Which state has more energy? It is reasonable to guess that state B needs more en- ergy. To prove this and study the distribution of the hard device energy function along the lines of constant 0, we take partial derivative with respect to ya in (24) for constant0. If 0 is a constant, then 7L = 1(0) and 7,, = F,,(0) are constants too. Therefore, A12 = A;2(y0) and s s’(y0) are also functions of ya only, and it is found that 8A12 9?. = 412m = (www—Wrrrmwns‘tv.) (55> O From (25), a similar process gives as t 1 __ = s ' = . (56) 870 a (70) F,(0) - 1‘”(0) Finally, substitution from (56) into (55) yields BA12 . , WHO" 11(0)) - W,(l" 1(0)) 37; a = A12 (yo) = I‘”(0) — F,(0) (57) This means that the hard device energy function A12 monotonically increases along each 39 line of constant 0 as yo increases. This also implies that (0, yo) = (0, I‘l(0)) gives the minimum value of A12 for a given constant 0 whenever 06 < 0 < 1. In conclusion, state B indeed has more energy than state A. It will be useful in what follows to define the ratio appearing in (57) as a pseudo stress function _ WIKFIA‘S» " WAF 1(0)) Ops ‘ r,,(o)-r,(o) a 0108(0), 0c < 0 S 1. (58) Equation (58) means that 0 p s is the average height of the stress-strain curve in the interval 8A y1 = 1‘,(o) 0 C and 0 = 1, however the limiting values of 0ps(0) can be ob- tained by the following: Taking the limit 0 —) 1 from (58) gives . _ . W170" ”(OD-WAY 1(0)) J‘fl‘f’ps‘“) ’ ell—Tr r,,(o)—r,(o) ’ (59) It is to be noted that W,’(1‘,(o)) = 017(0). W,,’(1“,,(6)) = ol“,,’(o). (60) Using the L’Hospital’s rule (since now (58) is of the form %) and substitution from (60) into (59) gives 1' 0 ( )— 1' or”’(6)_cr”(0) — 11m0 — 1 (61) Ol—rinl P50 — 01-311 l‘,,’(0)-l‘,’(0) - 0—)l _ ° Similarly taking the limit 0 -—> 0C from (58) and using the L’Hospital’s rule (since now (58) is of the form 2) gives 4o . _ . W,,0, floc<0<1, ( (63) 0ps(0)=0, if0=lorif0=0c. An example of the function 0ps(0) is shown in Figure 13. Moreover the first derivative of 0 p 3(0) is given by 0' ’(O') _ (Wll(rll(o)) — WKFKO’» " O (Fl/(O) " FAG» ) (IE/(0') - rII’(G)) ’” (rubs—13(0))2 = ( 6mm) - 0) (17(6) - I‘,,’(o)) r,,(o) — r [(0) (64) Note from (63)], (39)2, and (43)2 that 0 P 3(0) > 0, 17(0) = 1 and I‘ll’(0) < 0 in the inter- L 0.5: cc 0.6 0.7 0.8 0.9 O‘ 1 0.5 Figure 13. The function 0P,( 0) with the choice of the function 0,,(7) given by (a.1) of Appendix-A. Note that (63), ensures that this curve is above the dotted line. 41 val 0c<0< 1. It is also to be noted that 1“”(0)>1‘I(0) in this interval. Therefore 0ps’(0) > 0 in the interval 06 < 0 < 1. To complete this subsection we now look at the two limiting cases: 0ps’(0c) and 0ps’(l). Taking the limit 0 -> 1 from (64) and using the . L’Hospital’s rule (since now (64) is of the form %) gives l (0,40) - 0) (17(0) - F,,’(o)) 011-11161" (0'): 0111-111 I‘,,(0) - 11(0) (Op,(0) - 6) 11"”(0) —- 1‘,(0)' Note that there is a finite value in the first term of (65), while the second term of (65) is of (65) = 0111111 (F, (G) - r” (0.))011-1‘;l the form 9 . Therefore, using the L’Hospital’s rule gives 0 lim 0 (0) lim (1‘ ’(0)— (0)) lim (0 ”‘1 o)- I) = r” 11"” (0')- F, (0') (66) = 1- lim 0p '.(0) Therefore 1im10ps ’(0) = - (67) Similarly from (64) the other limiting case of 0 —> 0 c can be shown to be given by h , (0,40) - 0) (17(0) - 1‘,,’(0)) 0300’" (O)=0—>0 FII(O)-FI(O) , (68) . I“, (0) " I‘11 (0) lim = hm 111110(O'p (0)- 0). 0 40c 1.1/(O)— [(0’) 0—> Note that F,,(0) — 1“,(0) is divergent as 0 —> 06. Hence let’s assume 1‘”(0) - 17(0) has an n-th order pole as 0 —) 0 c. This implies that lim (I‘ll(0) — F,(0)) (0 — 06)" = constant (69) O -) 0‘ . and FH’(0) - 17(0) has an (n+l)-th order pole: 42 lim (r,,'(o)—r,’(o)) (o-oC‘)”1 = (-n) constant. 0—>0‘ Therefore 1‘,’(0)- Fifi“) _ . n 0—90c Fl,(0) - I‘,(o) _ 0—>0CO'- O’c. Substitution from (71) into (68) gives (0ps(0) - 0) n G- 0c lim 0 '(0) = lim 0-—>0c P5 o—io Note that (72) has the form 3 and using the L’Hospital’s rule gives that 1' ’ = ' ’ - . 013100” (0) Olgnocmps (0) 1)n From (73) and 0ps’(0) > 0 whenever 06 < 0 < 1 it is found that n-l 00 if n51, lim 0 ’(0) = 0-90‘ P5 { n if n>l, where n implies the rate of 0"”(7) —> 0c as 0 -> 06 by virtue of (69). (70) (71) (72) (73) (74) If the function 0,,(7) given by (al) of Appendix-A is considered, then it is found that n = l and lim 0ps’(0) = oo. Similarly, if the function 011(7) given by (b.1) of Sec- o—ioc tion Appendix-B is considered, then it is found that n = 1/2 and lim 0pS’(0) = 00. In 04¢ summary ( 0ps’(0)>0, if 06<0< 1, 0p5’(0) = 0.5, if 0 = 1, n Ops’ky) = {n _ 1 oo ifO=OcandO1, k (75) The first derivative of the function 0ps(0) with the choice of the function 011(7) given by (al) of Appendix-A is shown in Figure 14. 43 6.3 The hard device energy function for 912 on the lines of constant 7,, From Section 6.2 we found that the hard device energy function A12 monotonically in- creases along each line of constant 0 as 70 increases. Now we study whether the hard de- vice energy function A12 along the lines of constant 7 0 is increasing or decreasing with 0. Consider two different states with the same strain but different stresses on any line of con- stant 70: state A corresponds to (0, 70) = (0 A, 7A); state B corresponds to (0, 70) = (03, 73) = (08, 7,) where 03 > “71' Is it true that state B has a larger energy value? To answer this question and study the distribution of the hard device energy function along the lines of constant 70 again, we take partial derivative with respect to 0 in (24) with fixed 7 a. If 7 o is a constant, then A12 = A'l';(0) and s = s" (0) are functions of 0 only. From (24) and (60) it is found that 3A12 93 7. = A370) = (w,(1",(o)) - WrAr11(°)))sm(°) + (76) 0l‘l’(0)s + 01‘,,’(0) (1 - s) . Similarly, by the same mathematical process in (25), it is found that I O L r 1 r 0.5 = CC 0.6 0.7 0.3 0.9 0' 1 Figure 14. The first derivative of function 0p,(0) with the choice of the function 0,,(7) given by (a. 1) of Appendix-A. as n , I‘,’(0)s + I‘,,’(0) (1 — s) so (°’ = (r (o) - r (0)) (77) 7., II I Substitution of (58) and (77) into (76) shows that aAIZ 9‘ I I I 50 Y = A12 (0) = (1‘,(0)s+ 1‘” (0) (1 -s)) (0- 0ps(0)). (78) According to (34), equation (78) is equivalent to 3A” - " (o — o (0)) (79) 56 'Yo - 56 5 PS . 3A12 Since (63)] gives that 0 - 0 (0) < 0 whenever 0 at 0 or 0 at 1 it follows that and pr c 53 1 ——0 have opposite signs for 06 < 0 < 1 . This also implies that any extremum solutions of 8 80 the hard device energy Alzalong the lines of constant 70 occurs where 7 0 is stationary (stops increasing or decreasing) on the curves of constant s as a function of 0. Recall from Section 5 that (0, 70) on $2 is a minimum point in the (0, 70) -plane for the correspond- ing curve of constant s. Therefore it can be concluded that points on Q2 are also the ex- tremum points of the hard device energy function A12 along the lines of constant 7 0. The hard device energy function A12 at (0, 7 a) on E2 could be a maximum or a minimum along the lines of constant 7 0. This implies that , , o aAl2 rl(0)s+r”(0) (l-S)=O¢9§Es=0¢=956 7:0. (80) Equations (77) and (78) show the following relation 3A a 3;” = 5% (o-o,,(o)) (rum) —r,(o)). , (81) Y, Y. Similarly since 1‘”(0)— 1‘,(0) > 0 and (63)1 gives that 0 - 0 p 3(0) < 0 whenever 0 at 0 c or 45 3A 0 at 1 it follows that — 80 12 and 3%, have opposite signs for 0 c < 0 < 1. Although equa- ‘Y, 7 tion (79) and equation (81) are similar, they have different meanings. Equation (81) implies that any extremum values of the hard device energy function A12 along the lines of constant 70 takes place where s is stationary on the curves of constant 70 as a function of 0. To de- termine whether these extrema of energy are a maximum or a minimum we calculate the second partial derivative of A12 with respect to 0 from (78) and then using (80) gives 2 3 A12 302 9‘ II 2 = A12 (0) = 133770 (cs-omen. (82> 7. Again since (63)] gives that 0-0ps(0) <0 whenever 0¢0c or 0¢l it follows that 32111 802 2 82 311d 5370 7. 5 have opposite signs for 0c < 0 < 1. Since, from Section 6.2, 70 has a minimum on the lines of constant s at points on Q2, it follows that A12 has maximum on the lines of constant 70 at $2. Thus $2 is not only the set of ordered pairs (0, 70) which are minimum points to the family of constant 3 curves, but 32 also gives the set of ordered pairs (0, 70) which are maximum points to the energy A12 along the family of constant 70 lines. Returning to the question of which state A or B on the line of constant 7 0 has more energy we see that if 70 2 1 then state B has more energy than state A since the hard device energy function A12 increases monotonically with 0 whenever 70 2 1. But if 70 < 1 then since the hard device energy function A12 will increase monotonically to a maximum point on $2 and then will decrease monotonically. It follows that state B is not necessarily car- rying more energy than state A. Let’s define Q3 to be the set of ordered pairs (0, 70) with 0c < 0 s 1 that are not 46 pure phase-I solutions but which have the same energy A12 as the pure phase-1 solution at the same constant 7 0. Then it is found that $3 = {(0, 70)[70= (2°ps(°)'°)’°e<°5 1}, (83) where 0 p 3(0) is defined in (58). Therefore, at a given constant 70 obeying 0 c < 70 < 1, the hard device energy function A12 at points on El and G3 are equal, and this energy is less than the hard device energy function A12 at points on $2 (see Figure 15). With the help of (53) it can be concluded, at a given constant 70 obeying 0c < 70 < 1, that < A12] - A12] < A - . 84 on $1 on $3 12Ion $2 ( ) 12'0“ O = O" 6.4 The hard device energy function on the curves of constants We now look at he behavior of the hard device energy function A12 on the curves of con- stant s. Recall from Section 5 that there are two types of behavior to the constant s curves in the (0, 7 a) -plane: monotonically decreasing (s S s ) and decreasing then increasing trans (s > s ). For the first type it is found that A12 also decreases on the constant s curves as trans 0: 0c —> 1. In the second type with the help of (a.14) or (b.14) it is found that A12 also de- creases on the constant s curves as 0: 0c —9 00 (s) but then A12 increases on the constant 2 s curves as 0: 0¢2(s) —> 1 . To show this we take the partial derivative with respect to 0 in (24) with fixed s. It is found that A12 = A;;(o) and from (24) and (60) 8A 80 12 = A370) = (I‘I’(0)s + r,,'(o) (1 -s) )0. (85) S From (85), (34), (79) and (81) we find the following relations. —&YO s—fis 3A 0=E 8A 3'6 12 12 0’ _ as W _ 33700(F”(0)-FI(0)). ‘ (86) Ya 47 AI ooooo collar:III'IIOIIOIIIIIOOA yo = I“11(6) IIIIO'IIIIII'OO'IIIIIIII vIIIIIIIIIIIIIIIIIIIIIIIIII’IIIIIIIIIIIIIO’IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII I'O’O‘IIO'OO'IIOIO'II'II'Ili'IA / / J / /. __ /,../.. v... I t IIIC'IfIvl'IOOIII': I'I'IIIOI s AAtl l 8". A IIIIII"IIIIIL ['0’ A t AA 7 A-A profile 2 1 ¢ Qcl; \‘Illololollllll II”” lllll ’l’l’ IlllIOOIIIIaIIIA lllll 0' e B-B profile 70 < 1 C Figure 15. The hard device energy function A12 in the solution region 912. 48 Since 0> 0, 0 - 0ps(0) < 0 and 1‘”(0) — 1‘,(0) > O in the interval 06 < 0 < 1 it follows BA 81‘” and as have the same 51 ns while . {To 7 g a; that 56 0 ’3—0 3 12 has the opposite sign. This 7 0 implies fi'om (85) and (86) that , . 8A12 310’ 8A12 as Meanwhile, from (85) we can obtain the second partial derivative of A12 with respect to 0 and then using (87) gives 82A 2 12 _ tee” _ a “—802 - A12 (0) -— —aoz70 0. (38) S From (88), (82), and (77) it is found that 3271 2 32A 32 12 a 12 O s =— O'=-— —=—0I‘0—I‘0. 89 802 3 ad 70: 802 o_o-ps(o) a 2 ( ”( ) [( )) ( ) It is important to note that 0 > O, 0 — 0p5(0) < O and 1‘”(0) - 1’,(0) > 0. Recall from Sec- tion 6.3 and with the help of (89) that A12 and 70 on the lines of constant s are minimum at $2, and that s on the lines of constant 70 is minimum at $2, and that A12 on the lines of constant 70 is maximum at $2. Therefore $2 is an interesting curve. It not only gives the location of the minimum energy and the minimum 7 0 value on the curves of constant 3, but it also gives the location of the maximum energy and the minimum s value on the lines of constant 7 0. 6.5 The minimum energy criterion for the hard device Recall from Section 6.3 and Figure 15 that if 70 2 1 then the hard device energy function A12 monotonically increases as 0 goes from 0 c to 011(70). However if 70 < 1, then the 49 hard device energy function A12 will increase monotonically to a maximum point as 0 goes from 0 c to 002(70), but then it will decrease monotonically as 0 goes from 0% (70) 2 to 70. The minimum energy solution for the hard device will be the 7 a = 1‘[(0) = 0 line for 0 < 70 S 0 6. However according to the profile’s shown in Figure 15, the minimum en- ergy solution would be the 0 = 0C line for 70 > 0c if this line was attainable (see Figure 16). In fact for a given 70 there are 3 situations: 0 < 70 S 0 c, 0 c < 70 S l and 70 > 1. In the first situation, 0 < 70 S 06, the solution is unique. In the second situation, 0c < 70 S 1, every point on the line segment 06 < 0 < 61(70) of the constant 70 line correspond to a two-phase solution, while 0 = 7 0 is a pure phase-I solution. In the third situation, 7 a > 1, every point on the line segment 06 < 0 < 011(70) of the constant 7 0 lines also corresponds to a two-phase solution while 0 = 0" (70) is a pure phase-11 solution. For the second and third types there are an infinite number of possible solutions even for the (1,2)—case, however the values of A12 varies among these candidate solutions. Under a minimum energy solution criterion if the bar is elongated at its end with a quasi- e \“ .e . .\ h e e .n e «/ rim) ' V 0' 1 O’ Figure 16. The minimum energy solution for the hard device. 50 static monotonically increasing displacement 50 = 50(t), 50(t) 2 0 (a more detailed dis- cussion of quasi-static loading is given in Section 8), then the loading diagram will be as shown in Figure 17. This figure shows that under a minimum energy solution criterion the material first behaves as if it were linearly elastic and then behaves as if it were perfectly plastic for the hard device loading. However, unlike conventional plasticity, unloading af- ter 70 > 0C (corresponding to yielding in plasticity) simply retraces the path in Figure 17. ‘0 o o I o '0 .1 'o O I ‘0 0.. o..- ..... 0" '9 ................ 1 7, Figure 17. The loading-deformation diagram for the hard device under a minimum energy solution criterion. 7. The soft device energy function on (212 Section 6 was concerned with the analysis of the hard device energy function Anon the solution region. In this section we conduct a similar analysis for the soft device energy function 312. According to (26), the soft device energy function 512 is given by -— A u- —— e—I —-—a 12 1210’ Yo) = W,(F,(0))312(0, 70) + W”(1‘”(0)) (1 - 312(0, 70)) — 070. The corre- sponding level curves of fixed 312 in the region (212 are also functions of 70 and 0. An example of these level curves is given in Figure 18 for the same particular material used for the constant 5 curves in Figure 10 and for the level curves of fixed A12 in Figure 12. For this material and each fixed value of 512 e (—0.5, oo) , the equation sum, 7 a) = 512 de- fines a smooth curve on the (0, 70) -plane. There is a family of such curves in the (0, 7o) - plane which correspond to the different values of 312 e (-O.5, oo) . This section is orga- nized similar to that of Section 6, we examine the soft device energy function 312(0, 70) on the boundary of 912, on the lines of constant 0, on the lines of constant 70 and on the lines of constant s. 7.1 The soft device energy function on the boundary of On The soft device energy function 312 on the boundary of (212 is given by 02 = -7; on 70 = 1‘,(0), '2" 2’ - i W”(r 11(0)) " 01"”(0), on 70 = I31(6)’ (90) - -O' 0'2 70 c ’ LWKFKOCD + W,,(1‘”(0c)) “1(0) :0 - 0670 = -?, on 0 = 06. C i W,(r,(o)) —02 = — [11 l 12 The boundary 7 0 = 171(0) of the solution region 912 gives the pure phase-I solutions. As a result, the soft device energy function is 312 = -02 / 2 . The boundary 7 0 = 1"”(0) of the 51 52 3 l l V :12 = 0.06 70 312 = 0.02 512 = —0.02 a], - -0.06 512 = -0.10 512 = -0.14 2 .— 312 = —0.18 512 = —0.22 512 = —0.26 on = -0.3 -0.34 -0.38 —0.42 1 i = _/ -0.46 /‘ -12 _ —0.125 ‘\ an = —0.5 512 = —0.125 0 J I 0 0.25 0.5 = 0' 0.75 0' 1 C Figure 18. The level curves of constant soft device engrgy function En in the solution region 912 with the choice of the function 0,,(7) by (a.1) of Appendix-A. 53 solution region 912 gives pure phase-II solutions. Consequently, the soft device energy function is 312 = W,,(l‘”(0')) — OI‘”(0') and it will depend on the choice of the function 6,,(7). We now show that the soft device energy function E. on 0‘ = O'c, which is the 12 non-attainable limit, is given by 312 = -c:/ 2. From (26), (25) and (6) the soft. device en- ergy function 512 is given by _. 7__L :12 = __TL- -070 on o = 0c. (91) Y =°c 711"” With the help of (51), (52) and (91) it is found that ._ 0. 63 :12 = CC (YO-7) _OYO = -_2—° (92) It is to be noted that if y = 0' , then E. = -O’2/2 which indeed matches the energy for o c 12 c the pure phase-I solution. 7.2 The soft device energy function for Q12 on the lines of constant 0 Having investigated the soft device energy function 512 on the boundary of solution region Q naturally we now study the soft device energy function 512 inside the solution region 12’ (212. We first explore the lines of constant 6. Consider two different states with the same stress but different strains: state A corresponds to (o, 70) = (0A, 7A) ; state B corresponds to (0', yo) = (GB, 73) = (0A, YB) with 73 > 7". To determine which state has a larger en- ergy value we take partial derivative of 512 with respect to ya in (26) with constant 0'. Since 0 is a constant, 7L = F,(o) and 7R = I‘”(o) are constants too. Therefore, 12 -'-‘-12(Y Yo) and s = s (70 ) are also functions of 70 only, and it is found that - 54 8:1 570 2 = sigma) = (W,(r,(o» — W,,(r,,(o») s‘m) - o. (93) 0 Finally, substitution from (56) into (93) yields 8312 WHO" 11(0)) ‘ W,(I‘ 1(0)) 370 o _ mo) - mo) " ° = 0113(6)- 0 > 0' (94) Equation (94) shows that the soft device energy function 512 monotonically increases with 70 along the lines of constant 0' for cc < o < 1 so that the value 70 = I‘,(0') always gives the minimum energy solution for a given constant a. In conclusion, state B has more ener- gy than state A. It is to be noted from (94) and (57) that the energy functions A12 and {-312 of both devices monotonically increase with ya along the lines of constant a for ac < o < l. 7 .3 The soft device energy function for 912 on the lines of constant 70 We now explore the lines of constant 70 and investigate whether the soft device energy function 312 on the lines of constant 70 is increasing or decreasing with 0. Recall from Section 6.3 that we considered two different states with the same strain but different stress- es: state A corresponding to (6, yo) = (0A, 7A) and state B corresponding to (0', yo) = (08, YB) = (03, 7A) where OB > (3". Do we find a conclusion similar to (84)? To answer this question we take partial derivative with respect to o in (26) with fixed 70. For a given constant 70, 512 = 33(0) and s = s" (0') are functions of a only. From (26) and (60) it is found that ._ =.=..“'(o)= (w (cm—W (o»)s“'(o)+ ac 1° 12 ’(r' ”W” . (95) or,'(o)s" (o) + or,,'(o) (1 - s“ (0)) — yo. 55 Equations (58), (95) and (77) show that 30 2 = (I‘,'(o)s”(o) + r,,’(o) (1 — s" (on ) (o — 0,10» - 1,. (96) 0 According to (34) and (86) equation (96) is equivalent to 3:512 86 = - 7,, + E0 (o - 0,40» a" ‘ (97) ‘ (r,,(o) — 13(0)) (0 — opp» . Y. = —yo+—a—6 3A Recall from (79) that in the hard device 5512 day” 3“ a? Y. 3 3:2 . 12 < . 3—0 0’ 6mm) However now for the soft dev1ce have opposite signs provided that and 0 might have opposite signs 70 36 3 or might have same sign. For example, since 0 < Ups and 70 > 0 if we have 330 > 0 then we find 6312 < 0 But we find 3312 36 7 ° 33 < 0. Also recall from Section 5 3y >0 onlyifa—é’ 3 7. that since the distributions of the lines of constant 3 are dependent on the material only as shown in Figure 10, the distributions are typical for both devices. Therefore the points on E2 give the minimum values of s in the solution region 012 on the lines of constant 70 which are the same for both the hard device and the soft device. Also recall from Section 5 and Section 6 that the points on $2 give the maximum values of hard device energy func- tion 312 in the solution region (212 on the lines of constant 70. But now the points on 32 do not give the extremum points of the soft device energy function 312 on the lines of con- stant 70. From (97) it is found that < 0. i (98) 56 This implies that the possible extrema of the soft device energy function 312 along constant 70 lines can only occur where the curves of constant 5 have negative slope in the (0', yo) - plane. Hence it is to the left of the curve 92. Similar developments can be found in the re- BE. lation between —— as 86 and —— 7. Be 12 7. To determine whether any such extremum of the soft device energy function 512 is a maximum or a minimum we calculate the second partial derivative of the soft device en- ergy function 312 with respect to 0. Thus (97), (98) and (64) give 32312 .... 802 = :1 (o) = 7. (99) 323 r as F 1' $3 (I‘ll(o)- ,(o)) (o-ops(o)) +337 ( ,,(o)- 1(0)). 7. ° Substitution from (89) and (86) into (99) gives 2: a -12 802 a 2 0 - ~12 (0') = 56—270 (CS-0,140» +3—O' r, s s — _ <0, (100) a2 87,, where from (46), (63)1 and (98) we find a 270 > 0, O’ < ops(o) and 33 0' < 0. Therefore S 83 if the point (0, yo) satisfies 37312 = 0 then the point (0, 70) gives a maximum soft de- 7. vice energy function value 512. It is easy to show from (90), at a given constant 70 obeying ac < 70 s 1, that the en- ergy on o = 0c is greater than the energy on a = 6,00. For example the difference be- tween the energy on o = 0c and the energy on 0' = 51(70) is given by 2 2 cc)_[_Y:) = “59.350, (101) 12lono=5,(70) — ( 2 2 — — u- _— h—t h—o 12'0“ O = 0c 57 Comparing this result to (53) it is seen that the soft device behaves in a completely different manner. As for a comparison of the energy :12 on O’ = 0c and on o = 011(70) for 70 > 1, examination of Figure 19 shows that neither a = cc nor 6 = 61,00) is ensured to have less energy on the lines of constant 70 > 1. Further investigation of this interesting issue (which is a possible topic for future study) might yield rules for constructing profile curves similar to those in Figure 15 but for the soft device. 7.4 The soft device energy function for (212 on the curves of constant s Recall from Section 5 that there are two types of behavior to the constant 5 curves: one has 70 monotonically decreasing as a function of O (s S s ) and the other has 70 decreasing trans as a function of 0 then increasing as a function of 0' (s > strans ). It is found for both types that the soft device energy function 312 decreases as a: 0'6 —) 1. To show this we take the _tt partial derivative with respect to o in (26) with fixed s, and it is found that 312 = :12(0') and 6312 _tt, i 35 = :12 (o) = —yo<0. (102) 3 Yo I‘”(o) V CC 16 Figure 19. The minimum soft device energy solution on the constant 7,, lines for ya > 1. 58 7.5 The minimum energy criterion for the soft device Recall from Section 7.2 that the soft device energy function 312 monotonically increases with 70 on the lines of constant 0. Therefore the minimum energy solution for the soft de- vice will be the line 70 = l“,(o) which corresponds to a pure phase-I solution (see Figure 20). Therefore under a minimum energy solution criterion if the bar is elongated at the end with a quasi-static monotonically increasing displacement o = 6(t) , C(t) 2 0, then the loading diagram will be as shown in Figure 21. This figure shows that the material behaves as if it were linearly elastic up to the maximum load value a = l. Figure 21. The loading diagram for the soft device under the minimum energy criterion. 8. Quasi-static motions with a single phase boundary 8.1 Definition We now turn our attention to quasi-static motions of the bar. Following Knowles [1979K] and Pence [1991PP] a quasi-static solution to either the soft device problem or the hard de- vice problem is a sequence of equilibrium states which at each instant t obeys the boundary conditions. Note that no inertia effects are considered here. Now the time t is being intro- duced and it plays only the role of a history parameter. Thus in the hard device 700) is pre- scribed and is assumed to be a continuous function. At each instant t if 0 < 70(t) S 0c then there is an unique equilibrium solution a = 70(1). Moreover at each instant t if 70(t) > Cc then for both the (12)-case and also for the (2,1)-case there is a family of equilibrium so- lutions parametrized by the instantaneous value of spatially constant stress 0’ (see Section 5). The parameter 0 will thus vary with time so that a = 0(t) and 6(‘y(X, t)) = 0(t), where 1(X, t) is the associated strain field. If we only consider either the (l,2)-case or the (2,1)- case then once the (unprescribed) 0(1) is found, it will determine a particular quasi-static solution at each instant I. According to (23) it is necessary that 0(t) be continuous for the location of the phase boundary s(t) = 311(00), 70(1)) to vary continuously with respect to time I. If we imagine a bar is in pure phase-I in the beginning and that later a phase boundary is created at X = L which then moves from the right end to the left end, then we have the following different cases during the motion. At first the bar is in the (1,1)- case; and then the phase boundary is created at the right end and moves to the left mean- ing that the bar has changed from the (1 ,l)-case to the (12)-case; then if the phase bound— ary hits the left end then the bar changes from the (l,2)—case to the (2,2)-case. From the above observation we can find (1,1), (1,2), (22)-cases during the motion. This set of cas- 59 60 es which gives the motion of a single phase boundary is the combination of the (1,1), (1,2), (2,2)-cases. The other set of cases which also gives the motion of a single phase boundary is the combination of the (1,1), (2,1), (2,2)-cases. Here the phase boundary is created at the right end and, at least initially, moves to the left. Note that there can be no direct switching between the (l,2)-case and the (2,1)-case. The switching between those two combinations is only through an intermediate pure phase state, either the (l,l)-case or the (2,2)-case. Similarly, in the soft device, 0(t) is prescribed and is assumed to be a continuous function. At each instant t whenever O < 0(t) 5 ac then there is an unique equilibrium so- lution yo = 0(t). Moreover at each instant I whenever cc < O’(t) < 1 then there is a family of equilibrium solutions parametrized by the instantaneous value of elongation 60(t) = 70(t)L at the end X = L. The parameter 70 will thus vary with time so that yo = 70(t). Once the (unprescribed) 70(t) is found it will then determine a particular quasi- static solution at each instant t. It is also required that 70(t) be continuous for the location of the phase boundary to vary continuously with respect to time I. Consider that the bar is elongated with a quasi-static load 70(t) starting from (0, yo) =(0,0). For the combinations of the (1,1), (1,2), (22)-cases of solutions, since the solution at each instant is represented by a point in the solution region (212, a continuous quasi-static motion is represented by a continuous curve in this region. For a prescribed monotonically increasing 70(1), any continuous curve with increasing 70(t) in this region represents a possible quasi-static motion. Therefore we can find many possible curves for the motion. Typical curves for such motion which are given as path (a)-(f) in Figure 22. In what follows we explain the meanings of these paths: I (a) The solution path will be restricted to 70 = F,(o) = o of the pure phase-I solutions so long as yo(t) < cc. 61 (b) The solution path could stay on yo = I‘,(o) = o of the pure phase-I solutions after 70(t) > Cc until 70 = l. (c) Alternatively, after 70(t) > cc but before 70(t) = 1, the solution paths could move into the interior of the solution region by creating a phase boundary at the right end. ((1) There is another possible path which will first move along the line 70 = l‘,(o) = 0' of the pure phase-l solutions up to 70 = 0 = 1 and then could move into the interior of the solution region. Note that this point, (0', yo) = (1,1) is the only point in the solution region 012 which can be regarded both as a pure phase-I solution and as a pure phase-II solution. Upon leaving this point the strain field could have smooth transfer from this pure phase solution to a mixed state of phase-I and phase-H by cre— atin g an internal phase boundary. To see this start with the pure phase solution point 70 A 1“,,(O) (f) (e) (b) I‘,(o) (a) > 4 y 9c l c. Figure 22. Some possible quasi-static motions in the solution region (212. as 7,, increases from 70:0. 62 (0, yo) = (1, l) and imagine some fixed s = so between 0 and 1. According to Section 5 each curve of constant s in the solution region ends at (0', yo) = (l, 1) . It is to be noted that the strains on both sides of the internal phase boundary so are equal to l at the point (6, yo) = (l, 1) so that at first both sides of the internal phase boundary so have the common strain value of 1. Suppose that the strain 1R in x > so increases and at the same time that the strain 7L in x < so decreases in such a way that 511%) = 61(yL). This causes the solution path to move into the solution region 912 along the s = so curve. It was assumed in the above description that the internal phase boundary remains fixed at s = so, but it is also possible that the phase bound- ary could move. This occurs if the solution path (d) enters the solution region {212 by a path which does not coincide with the curves of constant s. (e) Another possible path stays on the curve 70 = I‘”(o). The solution will give pure phase-II solutions as long as the path follows the curve yo = I‘”(O). (f) The final possible path splits off from the solution path (e) at some yo > 1 and moves into the interior of the solution region 5212 by creating a phase boundary at the left end. 8.2 Dissipation and driving traction Following Knowles [1979K] the rate of dissipation D(t) in a quasi-static motion is the dif- ference between the rate of the external work E o(t) and the rate of increase of strain energy E o(t) , therefore: D(t) = 15.0) - 8.0). (103) As shown by Knowles [1979K], quasi-static problems with non-monotone stress response do not require that D(t) = 0 as it is in conventional linear elasticity. The external work is 63 Eon) = fio“’6(v)dv (104) and the rate of external work is 13.0) = 0700). (105) The stored strain energy is x 8,0) = jgwmdz = W010)» + W010» <1 - s) . (106) Note from (6) that 7L = [1(0) and 7R = 1‘1(0),wherei = I or II and j = I or [I . Hence the stored strain energy for all four cases could be obtained by equation (106). The rate of increase of the strain energy is then given by 5.0) = (W(vL(r))—W(v,.(t)))s+ (govt/01(0))” (-C%W(y,,(t))) (1 -s). (107) where s = gs; is the quasi-static phase boundary velocity. It is to be noted that o‘itwmm) = 611. -‘—’-W0 (0) = or . (108) dt R R From (107) and (108) it is found that 5.0) = (W(71(t)) — Wont)» s' + ohms + 01.0) (1 — s) . (109) Substitutions from (105) and (109) into (103) yield the rate of dissipation 00) = 07,0) - (W(71(t))- Wont») s' — (0110).) + 071.100 - s) ). (110) Taking the derivative with respect to time t in (23) gives 7'00) = (mos + You) (1 - s) ) + (710) - 11(0) s‘. (111) Therefore from (110) and (111) it is found that the rate of dissipation is given by 00) = s“ (W(7R(t)) - W010» - a (71.0) — 710)) ) . _ (112) It is to be emphasized that the rate of dissipation for all four cases as given in (21) could be obtained by (112), however D(t) = 0 for both the (1,1) and the (2,2) trivial cases since then 64 1,0) = 110). Letf = fij bedefined by f =f.-,- = W(I‘,(o)) - W0“ 1(0)) - 6 (13(0) - 13(0)) 11 = ( (WU-‘50)) - 013(0)) - (W(I‘i(o’)) _ 011(0)) ) . ( 3) Note that equation (113) can be rewritten as f=fij = [lpll = 171-1); (114) where [[0]] indicates the jump across the phase boundary at x = s(t) and p = W(7) - 01 which is the one dimensional case of the more general energy-momentum tensor discussed by Eshelby [1975135]. Substitution from (113) into (112) and with the help of (6) obtains the rate of dissipation D(t) = if. (115) Hence (103) can be rewritten as 13.0) + (-f) s' = 5.0) (116) which may be viewed as a work-energy identity. It means that the sum of rates at which work is being done on the bar by the external forces and done internally at the phase bound- ary x = s(t) , is balanced by the rate at which energy is being stored in the bar. Consequent- ly, -f may be treated as the traction applied by the phase boundary “on” the surrounding material, or similarly, + f can be viewed as a driving traction exerted on the phase boundary by the surrounding material (Abeyaratne and Jiang [1989A]). The driving traction is similar to the notion of the “force on a defect” first derived by Eshelby [1956B]. It is also similar to the force on the interface between two phases derived by Eshelby [197013] and discussed by Rice [197 SR]. It is easily shown that if there are numerous phase boundaries then (116) generalizes to Ea(t)+ 2 HM = Eta). (117) disooniiiruities From (58) it is found that (113) can be written as 65 f,-,- = (01,3(6) - 0) (5(6) - 11(0)). (118) Note for the (12)—case and the (2,1)-case that this gives f12 = (com) - 0) (r,,(o) - 17(0)) sine). (119) f2, = (0,40) - 6) (13(0) - 1“,,(o)) Ef21(0) = -ftz(o). In the (12)-case, (63)1 gives ops > o and I‘”(o) > F,(G) whenever cc < 0 < 1. Therefore it follows from (119)1 that the driving traction obeys f1:z > 0. However, in the (2,1)-case, it is found from (119)2 that f = f21 < 0. In summary f12 = -f21>0, if oo 0. (124) Recall form (39)2 and (43)2 that 170;) = 1 and r,,'(o) < 0. Therefore [32"(0) > 0 so that 1 0.5 = C 0.6 0.7 0.8 0.9 O. 1 Figure 24. The first derivative of the functions of f12(o) for the example function 511(7) given by (a.1) of Appendix-A for material A and by (b. l) of Appendix-B for material B. 67 [12(0) is a concave function. Since f21(o) = —f12(o) it follows that f21(o) is a convex func- tion. For the (1,1)-case and the (2,2)-case, both the strain 7L and the strain 7R are equal which gives f = f11 = f22 = O and hence the rate of increase of external work is balanced by the rate of increase of strain energy. This is to be expected, since for these cases an i11- temal phase boundary has no real physical significance. 8.3 Admissible direction of movement of the phase boundary According to the second law of thermodynamics under isothermal conditions, the rate of increase of strain energy E o(t) cannot exceed the rate of the external work E o(t) (Knowles [1979K]). Thus it follows from (115) and (119) that the dissipation rate D(t) = s‘fz O (125) which restricts the direction of the movement of the phase boundary. A quasi-static motion is admissible if the dissipation inequality (125) is obeyed during the motion. Hence it can be concluded that the directions in which the phase boundary can move are not unrestricted, but must obey unrestricted if f = 0, (126) 20 if f>0, 5) 50 if f<0. Consider the (l,2)-case, then from (120) and (126) it is required that s 20 if cc < o < l . However if the stress 0 = 1 then the driving traction f12 = 0 so that inequal- ity (126) does not restrict the sign of s at O' = 1 (see Figure 25). This restricts the possible paths as previously shown in Figure 22. In particular if the bar is elongated with monoton- ically increasing yo(t) then the path must follow the line s = l of pure phase-I solutions when 0 < 70(1) < 1. Thus path (c) in Figure 22 cannot occur. Thus the stress eventually reaches its maximum value c = 1. After this further elongation allow both pure phase-II 68 a ll __ O K 70 = r1183) = O D} II n—n N I N u—u l "\4 AIM r4 II #I U) 0 I 1' 1 o 0.25 0.5: O": 0.75 O’ 1 Figure 25. The admissible directions for the quasi-static _r_notions in the solution region {212. The constant s curves are for the function 011(7) given by (a.1) of Appen- dix-A, but the diagram is representative for any function obeying (5). 69 solutions (5 = O) and true (1,2)-type solutions. This is because the admissible directions on the curve 5 = O of pure phase-I] solutions allows the paths to split off and enter the in- terior of the solution region 012 (see Figure 25). The (2,1)-case is similar to the (l,2)-case, but now the phase boundary starts on the other side and moves in the other direction. Hence the admissible path diagram is essential- ly the same (Figure 26). From the above discussion it is to be noted that the solution paths will start from (0', yo) = (O, O) , follow the 70 = 0' line, and reach the peak (0', yo) = ( 1, 1) for both hard device and soft device loadings in the (0', yo) -plane. Moreover since there is no phase boundary before this peak is reached it follows from (116) that the rate of increase of external work is balanced by the rate of increase of strain energy. Thus, the path 70 = o of pure phase-I solutions is dissipation free. Note that for the same reason the path 70 = THUS) of pure phase-II solutions is also dissipation free. Finally we note for both the (12)-case and the (2,1)-case that the phase boundary moves so as to decrease the amount of phase-II and increase the amount of phase-l. This is related to the instability of descending portions of a stress-strain response function as found by other authors for different problems (see Abeyaratne and Knowles [1989A], Kikuchi and Triantafyllidis [1982K]). 8.4 Comments on the admissible paths In what follows, equation (126) will be required to hold at all times t. We first consider soft device paths in which 0(t) is not only specified but is specified to be constant for all future time t. We consider arbitrary initial values of 70 in the solution region (212 on the associ- ated line of constant 6. An interesting question is whether paths starting at this initial point are allowed to move to the boundary of the solution region 5212. For each constant value of 0’ there are two possible boundary points given by 70 = l",(o) or 70 = rum). Every point , 1 i l o 0.25 0.5: cc 0.75 O’ 1 Figure 26. The admissible directions for the quasi-static motions in the solution'region 9.2,. The constants curves are for the function 0,,(7) given by (a.1) of Appen- dix-A, but the diagram is representative for any function obeying (5). 71 V o 1 o C Figure 27. Admissible constant 0 paths that reach the boundary 70:1",(0') for the soft device. All points on the constant 0 lines can reach points on this boundary. in the solution region could move to the boundary point 70 = 1",(o) (see Figure 27) since a downward path is always admissible (see Figure 25). The points inside the solution region could not move to the boundary point 70 = F"(0') since the upward direction is not admis- sible (see Figure 25). Therefore the only point that can “move" to 70 = I‘"(o) on an ad- missible constant 0’ path is the boundary point 70 = I‘”(0') itself which is already there (see Figure 28). Similarly let’s consider hard device paths in which 50(t) is not only specified but is specified to be constant at all instants t so that 70 is fixed. The same interesting question is whether similar paths could move to the boundary of the solution region. For each constant value of ya, there are two possible boundary points given either by 0' = O", and ya = F,(o), or by 0' = 6c and ya = F"(o') depending on the constant value of 70. Those points inside the solution region with 6(70) > 691(7) and 70 s 1 could move to the bound- 72 cc 1 a Figure 28. Admissible constant 0 paths that reach the boundary 70=F,,(0) for both the hard and soft device. No points, other than the boundary points themselves, can reach the points on this boundary. ...... ‘ ..... I‘,(o) <— : E > V 0c 1 0 Figure 29. Admissible constant 7,, paths that reach the boundary 70:1“,(0) for the hard device. Those points which are on the constant 70 line with 5(7)>392(7) can reach this boundary. 73 ary point 70 = I‘l(c) (see Figure 29) since the rightward path is also admissible in this area (see Figure 25). Those points inside the solution region with 6(70) < 66570) and 70 S l, as well as those points inside the solution region with 70 > 1 could move toward the boundary point 6 = 0C (see Figure 30) since the leftward direction is admissible in these areas (see Figure 25). It is to be recalled that the boundary point 0 = cc is an unattainable limit as- sociated with minimum energy. Recall from Section 6 that the point (632(7), 70) on the curve Q2 represents an energy maximum on lines of constant 70. Thus admissible motions in Figure 29 and Figure 30 must move down the energy hills of A (see Figure 15). From the above discussions it is again seen in all cases that the final states on the boundary of (212 always like to reduce phase-H and extend the amount of phase-I. Because of the restriction in the direction of solution path we shall now briefly re- consider the six possible solution paths (a)-(f) previously depicted in Figure 22 and dis- cussed in Section 8.1. We again consider motions in which the bar is elongated with mono- tonically increasing applied displacement 70(t) > 0 starting from (a, 70) = (0, 0) (see Figure 31). (a) As before, the solution path is restricted to yo = F,(O) = U of the pure phase-I solu- tions until 70(t) > 0C. The bar behaves as if it were linearly elastic for both devices. (b) The solution path must now stay on yo = 13(6) = 0’ of the pure phase-I solutions after 70(t) > cc until 70 = 1. (c) After 70(t) > cc but before you) = l the possible solution path(c) of Section 8.1 can- not occur since this violates (125). (d) Other admissible paths enter the interior of solution region after the initial departure from (0’, 70) = (l, 1) . However the solution path cannot move arbitrarily, that is there is less variation since it cannot violate the directional restrictions in Figure 25 . E : J. v 0' l 0' C Figure 30. Admissible constant 7,, paths that reach the boundary 0: dc for the hard device. Those points which can reach this boundary are on the constant 70 line with oc< 5(yo)< Egg 1,) whenever oc<70<1, or with oc< 5(7)< 5,,(70) whenever yo 2 1. ‘ I I ' o 1 o C D Figure 31. Some admissible possible quasi-static motions in the solution region 912 as 70 increases from 70:0. 75 after it enters the interior of solution region. (e) The path stays on the curve 70 = I“, [(0) of the pure phase-II solutions. (1) The solution path enters the interior of the solution region (212. However the path must obey the directional restrictions in Figure 25 as do the paths (d). From the above discussions we find that the paths (a), (b) and (e) in Figure 31 are exactly the same as their corresponding paths in Figure 22 since there is no other choice. However the paths (c), (d) and (f) in Figure 31 are not necessarily equal to their corresponding paths in Figure 22 since those paths in Figure 22 do not need to satisfy (125) while these paths in Figure 31 must obey (125). Now that we have looked at the admissible paths, we return to the minimum energy solution diagrams in Sections 6 and 7. We find from (120) and ( 126) that, for a soft device case, the loading curve in Figure 2] is admissible. However, for a hard device case, the loading curve in Figure 17 is not admissible whenever O’c < 70 < l . This is because there is an energy hump on the curve Q2 (see Figure 15), so that attaining the absolute minimum on the line 0 = 0‘6 requires a climb over an energy barrier to leave the local minimum on the line 91. 9. The kinetic criterion Due to the non-uniqueness of the quasi-static motions solutions some additional criterion is needed which determines the unprescribed 0(t) for the hard device and which gives the unprescribed 70(t) for the soft device. Requirements of minimum energy provide one such criterion, but as mentioned at the end of the previous section, this can lead to violations of the admissibility condition (125). An example of a different type of criterion involves the use of a “kinetic criterion”. Abeyaratne and Knowles [1989A] have shown for problems with similar nonuniqueness that the quasi-static response of a prescribed force problem is then fully determined if one introduces the idea of an additional kinetic criterion. Here the kinetic criterion will be described by a relation f = (D(t). (127) where (p(s°) is a materially determined function that relates the driving traction f and phase boundary velocity S. Similar kinetic criteria were considered in [1988A], [1988AA], [1989A], [1990A], [1991A], [1991AA] and [1992A]. A necessary and sufficient condition to ensure that the dissipation rate D(t) obeys (125) is that q)(s) satisfies s' Cc. Our interest will be whether the paths eventually intersect the boundary of solution region (212 (recall Section 8) and if so, we also wish to know how long it takes and how fast the phase boundary is moving at this final time. For the soft device since the stress is given by o , S = kf12(o A) is constant throughout the motion and is greater than zero. Thus the phase boundary will eventually move all the way to X = L resulting in a pure phase-I solution (see Figure 27). To find how long it will take, the initial location of the phase boundary at point (0A, 7A) is needed. For the ( 1,2)-case, it is found that 1 3(0) = FII(OA)—YA _ rlI(UA)_‘YA. Won-170.)“ no.) “44’ (In contrast, if this was the (2,1)-case, then — -— 0‘ 3(0)_ 7A ,(0 A)) __ 7 ,( A), (145) r,,,,(0 )- r,(0A T(0 A) note that the sum of the two initial values of 3(0) from (144) and (145) is equal to 1). There- fore in the (1,2)-case from (144) the time t is given by : I -S.(0)_ I‘O.Al( O’-) 7A > (146) S kf12(0 A)T(0 A) (In contrast, if this was the (2,1)-case, then from (145) the time t is given by = '5 (,0) = W 0) 7" (147) S -kf21(oA)T(OA) Note that the formula for t is the same for both cases). Turning to the hard device there are two possibilities. Recall from Section 8.4 for each constant value of 7A , that there are two possible boundary points to the solution region 012. These are given byo = o ando = 7A whenevera c< 1. For the case in which either 6A < 602(7) with VA S l, or the case with VA > 1 (see the shaded region in Figure 30) it is not easy to solve this prob- lem as the velocity of the phase boundary increases through the motion. As an example for material A the velocity of the phase boundary is i = kf12(oc) = 00 when the solution path finally hits the line 0(t) = O'C while for material B the velocity of the phase boundary is equal to 5k/ 8 when the solution path finally hits the line 6(t) = Cc. If the stress is ob- tained from (143) then from (140) the time I when the phase boundary hits the end X = L is given by the following integral equation 1 5(0) + kjgf12(o(r))dt, for the (l,2)-case, _ (148) 0 = 5(0) +kjgf21(0(t))dt, for the (2,1)-case. Similarly for the case in which GA > 66 (VA) with 7A S 1 (see Figure 29), the velocity of 2 the phase boundary decreases through the motion and is given by i = kf12('yA) when the solution path eventually hits the line of phase-I solution at o = 6,01) = 7,1' The time needed is again given by (148). 9.3 Rate dependence of the admissible paths governed by the kinetic criterion In the previous section we briefly looked at the time needed for a (1,2)-case solution path to “drift” to the boundary of (212 when the prescribed loading condition (either 70 or o) is fixed. Now we turn to investigate the effect of different loading rates. We consider the sys- tem starting from the fixed point (0, 70) = (0, 0) in the solution region 012, and we try to find the kinetic solution for both the hard device and the soft device. For the soft device if 0 < 0(t) < 1, then the solution path will remain in pure phase- I (which is also the minimum energy solution, recall Section 7). While, for the hard device, if 0 < 70(t) < 1, then the solution path will also remain in pure phase-I (which does not give 82 the minimum energy solution for ac < o < I, recall Section 6). From the above observations the solution path cannot enter the interior of solution region (212 so long as either 0 < 0(t) < 1 or 0 < 70(t) < 1. Therefore let’s consider the sys- tem starting from the fixed point (a, 70) = (1, 1) at t=0 which represents a pure phase so- lution. Under further loading the solution path may now enter the interior of the solution region 912. We consider the hard device loading condition ‘th l+——=l+at for Ot, o where 7, and to are positive constants, and the loading rate a = 717/ to (this is shown as (c) in Figure 32). Thus (149) represents a family of loading cases parametrized by to and 7F. We are concerned with the effect of rate, that is the effect of changing to with 7,. fixed. If to increases, then the loading (c) is changed to some other loading, say (d); if t5 decreas- es, then the loading (c) is changed to some other loading, say (b). In particular, if to -) 0, then the loading (0) switches to the step loading (a). I claim that the solution paths will look like the six corresponding paths as shown in Figure 33. We now show this in the following. 1 + 7, ‘ A, . ”,.. /‘ 2 70(1) (20 ~““:(b) 37(0) ifs/(g), ,/ MAM"; ”I”, (c) f “,..-"v d, f «7” s .f s 1 ‘ .» (f) 4 5 > to (c) t Figure 32. The hard device loading condition (149) for difierent values to. Here 0 = to(a) < 10(b)<--- < to“) < to“) = 00 so that on = (10(3) > 0100)) >... > co“) > (100-) = O. 83 From (149) and for a loading with the loading rate a = 0 , the solution path follows the line of constant 70 = l as discussed in Figure 30. But for a loading with the loading rate a > O we have from (142), for 0 < t< to, that 7,0) = a = ('30) (1 + T'(G(t)) (1 - S(t))) - T(O(t))$(t) > 0. (150) From Section 5 if 70(t) 2 l and 6c < 0(t) < l we have that the slope of a curve of constant 3 in the solution region is negative (see Figure 9): 1+T’(o(t))(1—s(t)) = 8% <0 1170(021 and cc 0 and s’(t) > 0 whenever O", < 60‘) < 1, equation (150) now shows that c'r(t) < 0. Hence the slope of the solution path in the solution region (212 is negative whenever a > O and 70(t) 2 1. This means that the stress decreases while the strain increases for a hard de- vice loading with a constant loading rates a > 0. A I‘,,(o) l+yF ...................... . \ ‘ \ ‘ ‘ \ ‘ \ ‘ ‘ ‘ ‘ ‘ ‘ ‘ ‘ ‘ ‘ ‘ \ ~ \ . \ ‘ ‘ ‘ \ A ‘ > V 0 1 0 . C Figure 33. The relation between the solution paths corresponding to the loadings given in Figure 32 for the hard device. 84 Each loading (.. each (1) defines a different solution path 7(0) and to show this or dependence we shall put a in the argument 7(0, or). To show that a higher or path is above a lower a path (at each 0) it is enough to show that a“ > ab = 7(0, or“) > 1(0, orb) with 0 fixed. A formal argument establishing this ordering is given in Appendix-C, while here we give a shorter (but less formal) explanation. To argue that this is the case, note that we can- not have 7(6, (1"): 7(0, ab) at any fixed 0 <1 since this would imply that two different 7 curves pass through the same point. Thus either 1(0, (1") >‘7(O', (1”) or 7(6, a°)< )(O, ab) for all fixed 0 which a“ > (1". To find the correct direction to the inequality it is enough to consider a specific example, and so we take ab = 0 . However in this case 1(0, 0) =1 which is the lowest most possible path, so we conclude that 7(0, a“) >7(0, ab) when a“ > ab so that higher (1 curves are above lower (1 curves. We have argued that the rate of loading gives different solution paths. Thus differ- ent solutions are generated by different loading rates. Those curves would separate like in Figure 33. Note that some curves like (e) and (f) could hit 0' = 0c before t = to. Physically these results can be concluded by noting that there is a natural tendency for the solution path to drift down to the line a = ac, so that larger values of to allow more downward drift before 70 = l + 7],. For the loading (a), to = 0 and there is no time re- sponse for this drift Hence the solution path will then follow the path with the smallest pos- sible s which is the curve 70 = I‘”(o), 3:0, with l < 70 < 1 + 7F and will then follow the line 70 = 1+7], toward the line a = cc. For the loading (f), to = 00 so that there is enough time for complete drift response. Hence the solution path in this case will be the lowest most path which is the line 700) = 1. Finally the departure of the solution paths in Figure 33 from (0', 70) = (1,1) can be analyzed. This is done for Material A in Appendix-D. 10. Finite element analysis The theoretical analysis in Sections 2- 9 shows that the problem under study has nonunique solutions and that this nonuniqueness can be resolved by implementing an extra criterion, either one of minimum energy or else a kinetic relation. Now we study whether numerical analysis based on the finite element method can give the same results. Our problem is find- ing a state that satisfies the equilibrium equation (11) d = 152 _de (X) 0, 0 1 . Solutions in which all the elements are in phase-I are pure phase-I solutions. Solutions in which all the elements are in phase- II are pure phase-II solutions. All other solutions involve a combination of phase-I and phase-II elements. Nodes separating elements in different phase are phases boundaries. According to (14), the elastic secant modulus E function in the i-th element is given by 0(1‘) [(1') (i) __ E (i) . — (i) (170) “2 ’“1 Therefore the elastic secant modulus E (i) function in the i-th element is a constant. Equa- tion (170) can be rewritten as (i) (i) u -u ( 2 1 )< 51:1) [(1') - ’ (i) (i) (i) (i) E =E (u1 ,u2 )= i . . . (171) — 611((uéi) _u1( ))/l(‘)) (uél) _ul(l)) H — i ‘ - r . _ . (as) _ul“))/l(1) [(1) From (164) the element stiffness matrix K :10 becomes . (i) . (z) E l -1 (1) 1 -1 K . = _ = k 172 "’ 101-1 1] i—1 1] ( ) where the element stiffness km = —(7)— — . The force vector F ,f 0 is given by l “2(i)_u1(i) 89 - _ (i) Ff” = [0 ]. (173) 0(1) Finally (163) can be rewritten as (i) (i) _ E—(r[l ‘1‘“. = 0‘” [‘1]. (174) 1 ‘ -1 l u (1) l 2 Although the above equation (174), as well as many of the equations to follow, are written in a standard matrix form of a linear system, it is really important to remember that E‘0 in (174) is a function of up and uén so that (174) is a nonlinear equation. Therefore, the standard uniqueness and existence results for the linear theory might not hold for this equa- tion. In fact, since this equation of a discrete form comes from a continuous theory which has severe loss of uniqueness (before an additional criterion is considered), we would ex- pect that the discrete approximation provided by FEM would also have this property. 10.3 Assembly of element equations Since equation (174) is derived for an arbitrary element, it holds for each element from the finite element mesh. For an N «element case, we have the resulting system matrix equation [KJW+1)X(N+1)[U](N+1)XI = [Pi](N+1)x1° (175) Note that [K] is the global stiffness matrix and is given by k, -k, 0 0 0 0 -k1k1+k2 -k2 0 0 0 0 -k2 k2+k3... 0 0 0 [K] = (176) o 0 0 .../c,,,,2+k,,,_1 -k,,,_, 0 0 0 0 -k~_, kN_,+kN-k,,, o 0 0 0 -kN kN‘ 90 where the element stiffinesses k,- are given by k, = E (i) /l (i) . The nodal displacement vec- tor [U] is defined as T [U] = [“1 “2 “3 “Iv—1 “N uN+1] (177) while the nodal force vector is defined as T [F] = [6102 03017-1017 “rt/+1] 1 (178) where o, = a“) - (r(‘- + 1) is the force difference across the i-th node. Equation (178) is an assembly result of (165). In fact the forces on the internal nodes are equal to zero because there is no body force and the external forces are applied to the ends only. Therefore it fol- lows that am takes on the same value, say 0, in all elements and (178) is reduced to [F] = [—000 00 of. (179) From (169) and (177) the strain on the i-th element is constant and given by u. -u. y, = LF (130) l where the i-th element length 1,- = l“) . Since in this study all of the strains 7,. are positive, it is necessary that uj . (“HI-“1) 11' where l,- = l (i). Hence the stiffness matrix [K] and the nodal force vector [F] are each a function of the nodal displacements: 14,-. Thus equation (184) is a nonlinear system of equa- tions for the N -1 unknown values 112, ..., 11” which is simply written in the usual linear equation form. Note that the unknown value of stress can be obtained from (183) after solv- ing (184). For the rest of this thesis we shall limit attention to an FEM model in which all el- ements have equal lengths. Thus for an N -e1ement model the common length is l =L/ N . We now consider the possibility of solutions to (184) in which all the stiffnesses k,- are equal. Then by using Gaussian elimination one finds that the nodal displacement vector ‘ ... 92 must be given by “2 F 1 _ “3 2 5 = N' (186) uN_1 N-2 _uN, _~_1_ We then verify from (185) that (186) indeed allows for equality of all the k,- so that (186) is a solution to (184). Using (180) the strain 7, in each element is thus given by y, = 5/ (N1) = 5/L = 70. Here the strain in each element is the same constant so that all ele- ments are in the same phase. Thus (186) represents a pure phase solution to (184). Hence FEM always allows the possibility of a pure phase solution (we shall soon see that other solutions are also possible). This pure phase solution is a pure phase-I solution whenever 0 5 yo < l or equivalently 0 S 50 < L and is a pure phase-II solution whenever 70 > 1 or equivalently 50 > L. For the soft-device problem, the displacement at the first node and the stress in the bar are given by u 1 = 0, = a prescribed constant, (187) while the N nodal displacements: 112, ..., uN H are unknown. Similarly to obtain a new system equation we first insert the prescribed boundary conditions (187) into (175). After that we neglect the first row equation and move the ul term which is 0 to the right hand side of the equal sign. Finally the new system equation is given by "kl”, -k2 0 0 0‘ ' u2 ' 00- ... ... ... ... ... = ... . 188 O O ... kN-2+kN-1 —kN-1 O u~_l O ( ) O 0 ... 0 —kN kN uN+1 -O- ..n. I 93 Just as (184) was a nonlinear equation, so is (188). Now there are N unknowns 142, ..., 14,, + 1. Once again, a pure phase solution is available for (188), it is given by 112 P 1 ' “3 2 o = —. (189) “Iv-1 ‘ kl “N N-l _“~+1. _ N1 According to (180), the strain on each element is constant so that each element is wholly in one of two possible strain phases: either the lower strain phase (phase-I) or the high strain phase (phase-II). For a solution that is not pure, these phases however could be ordered in many different ways. For example, for N =4, these phases might be ordered as the following specific case [I 11 I I]. (190) on which the first element, the third element and the fourth element are in phase-I, and the second element is in phase-II. The number of phase permutations to (190) is 4 and the total number of phase arrangements is 24 = 16. As another example, for N =20, these phases might be ordered as the following specific case [Inrrrnnrrnrrrrrrurrunln]. (191) Here there are nine elements being in phase-I while there are eleven elements being in phase-II. The number of phase permutations to (191) is 21 and the total number of phase arrangements is 220 for a twenty-element model. In general the number of phase arrange- ments is 2” for an N -element model. This also can be clarified by Pascal’s triangle, which has the terms of polynomials of various degrees in the two variables I and II as shown in Figure 34. For a binary expression (1 + II)", we have 94 Number of Number of Number of Pascal’s triangle elements (N) cases (2”) terms (N+l) 1 1 11 1 2 2 I2 2111 112 2 4 3 13 31211 31112 113 3 8 4 I4 41311 612112 41113 114 4 16 5 (s=l) (s=.75) (s=.5) (s=.25) (s=0) Figure 34. The coefficients in Pascal’s triangle gives the number of phase permutations at fixed phase volume fraction for two-phase problems with different num- ber of elements. Here 2” gives the number of possible phase arrangements, and N+l gives the number of distinct volume fractions of phase-I to phase- II. N (1 +11)”: 2 C(N,m)1’"11”"", (192) 01—0 where N is the number of elements, m is the number of elements being in phase-I, N -m is l mix—7&7). is the number Of the number of elements being in phase-II and C(N, m) = possible permutations for which m elements are in phase-I and N -m elements are in phase- II. For equal element lengths and fixed m, all the C(N, m) permutations represent the same volume fraction of phase-I to phase-II. Note that pure phase solutions imply m = 0 or m = N in (192). Since C(N, 0)= C (N, N) =1 there is only one member in a pure phase permutation. For 0 < m < N the phase permutation corresponds to a multi-phase solution. Exactly one member of the permutation group has phase-I on left and phase-II on right and so corresponds to the (1,2)-case solution of Section 4. Moreover exactly one member of the permutation group has phase-I on right 95 and phase-II on left and so corresponds to the (2,1)-case solution of Section 4. The above two members are the only multiple-phase solutions with a single phase boundary. This phase boundary location is limited to nodal values and so must be at l 2 S = 0, l—v, l-V’ ..., l. (193) For the (12)-case solution the phase boundary location in (193) also gives the fraction of phase-I. Note also that all the remaining C(N, m)-2 permutations involve more than one phase boundary. It will be convenient to identify all members of the permutation group by its (12)-representative. Thus the specific case of (190) would be represented by the follow- ing general case (IIIHX (90 where <11111>=[11111],[11111].[11111].[111H]. (195) There are a total of N +1 possible permutation groups, denoted by < >. Each of these N +1 groups corresponds to a different phase boundary location (or equivalently volume fraction of phase-I) as given in the list (193). 10.5 Iterative solutions to the system equation To solve the nonlinear system equations (184) and (188) an iterative method may be cho- sen. Both equations can be rewritten as an iterative form in the following way [g(uj'_)ujm = 15(1),). ‘ (196) Note that the subscript j is the nodal number and the subscript i is the iteration number. The iteration procedure represented by (196) is the following: (a) Choose an initial guess (i =1) for the nodal displacements u j, (b) Compute the element stiffness fluj) using (185). Compute the nodal force vector 96 [(14)-i) from the nodal displacements uj‘ using (184) and (185) for the hard device problem or using (185) and (188) for the soft device problem. (c) Solve (196) for new nodal displacements u J. 1+1 ((1) Check the norm between 14,- 1 and u, by “ulna - “L“... = max (lulnt — “It" luznr — “21" H” lull/+1.41- uN+ Ill) ' (197) If the norm [1 u. < e, where e is sufficiently small then the problem is con- 1 ‘ 1.“... sidered solved and the process is stopped. Otherwise replace aim by "1'." go back to the step (b) and repeat again. In step (d) there are many possible methods to compute the norm. The best known of these is the Euclidean norm, such as used in the least-square method. But this is not our main con- cern, and we choose (197) to compute the norm. The method outlined above is the well known method of direct iteration or succes- sive approximation [19800]. Other more sophisticated iteration procedures are also possi- ble (such as Newton-Rhapson [19800]) which may have advantages from the point of view of convergence or convergence rate. This direct iteration procedure, however, will be suf- ficient for our purposes. From the above steps we expect to be able to find one solution from each initial guess (provided that the iteration procedure converges). But since the system equation, ei- ther (184) or (188), is nonlinear, it is possible that there may be multiple solutions. How can we find all of these solutions? This will be explored next. To find all of the solutions, the initial guess displacements uj are chosen from throughout the possible range. There are several restrictions to this range and thus to the initial guess displacements 14,-. For example all of the uj obey (181). Moreover for the hard device the displacement at N +1 node is equal to 5: uN H = 5. For the hard device, since 5 is given the possible range of the uj, could be from 0 97 to 5. But for the soft device, the stress 0 is given while 5 is unknown. However for a given a the strain 70 varies from F,(o) to 1“,,(0). Therefore the possible range of the uj, could be from F,(o)L to I‘”(G)L, where L is the total length of the bar. 11. Solutions of the system equation with N=2 11.1 Formulation System equations (184) and (188) for an N -element model were set up in Section 10. To better understand the relation between the solutions to those equations (possibly nonu- nique) and the initial guess displacements we will explore them for the simple case of 2 el- ements (N=2). For this simple case we also study how the energy functions depends on the nodal displacements. Under a hard device loading with 2 elements there are two prescribed conditions: it 1 = 0, 143 = 5, so that there are one unknown displacement uz and one un- known stress 0. Hence equation (184) is reduced to E E E _1, _3 [1),] = .35. (198) l l l where l is the length of each element and l = L/ 2. But under a soft device loading, the given conditions are u1 = 0, o = a prescribed constant, so that there are two unknown dis- placement uz, u3. Hence equation (188) becomes [El E2 E27 __+ _. __ u 1 1 1 2 = [0]. (199) T T] Both equations (198) and (199) are nonlinear since E1 and E2 are also functions of the un- known nodal displacements in the problem. For the hard device loading of equation (198), E1 and E2 are given by 98 99 E1 = E1042) = { 6,,(142/1) = __ > E” 142/l , 112/1- 1, E,=1, (5—u2)/l.<_1, E2 = E2012) = { 611((5-142) /l) E”: (5_u2)/l , (5-u2)/121. (200) For the soft device loading in equation (199), equation (200) continues to hold, provided that the prescribed value 5 is regarded as unknown. Thus one can replace 5 by u3 in (200) so that E2 = E2(u2, 143). For both loading devices it is easy to see that o = [2114/ I = E 2 (u3 - u2) /l . In what follows our main concern is the hard device loading, so that (198) holds. For each E ,- equation (200) gives two choices which are determined by the strains: yl=u2/l and 72: (5 - 142) /l. Therefore there are 22 = 4 possible forms to (198) corre- sponding to the phase arrangements [I I], [I II], [II I] and [II II]. The pure phase solutions (186) give 112 = 5/ 2 for the phase arrangements [1 I] and [II II]. In what follows the fol- lowing notations are the same: (1.1)= [1.1]. (1.2)= [1.11]. (2. 1)= [11.1]. (2,2)= [11.11]. (201) 11.2 Direct search for FEM solutions In this section we directly analyze each of four cases (201) individually. Such a direct search is feasible in this case (N =2) since these are only 4 possible phase arrangements to examine. In order to understand more about (198) it is convenient to define @042) to be the following “system function” so that solutions of (198) are given by roots of @012) = 0. If the following non-dimension- al parameters are introduced: 100 _ “2 5 u2 = r: 7, = Z. (203) then _ .. E270 9(112) — 142-m. (204) Applying the inequality conditions on the right hand side of (200) with 141 = O and 113 = 5 gives the following conditions on 112 and 5 for the four cases (201): (1,1)-case aquland 5—u2Sl=~5SL wyoS 1, _0 2’ 5 _ 70 (2,1)-casemu2>land5-u2Sl=>u2>-2- @u2>-i-, (2,2)-case¢=u2>land5-u2>l=>5>L «=70>1. (12)-case mquland 5-u2>l=>u2S§ ©172S (205) Note that 172 is restricted to 0 = 171 S 172 S 173 = 70 (see Figure 35). There are four differ- ent regions: region A for the (1,1).case, region B for the (1,2)-case, region C for the (2,1)- case and region D for the (2,2)-case. Figure 35 also verifies that both the (1,1)-case and the (2,2)-case cannot exist for the same 70. A D 70 B l A=(1,1) B=( 1,2) C=(2, 1) 1/2 =(272) Figure 35. The regions of four cases of possible 172 considering a given 7, for material A. 101 The roots of 9(172) = 0 in (204) are the solutions of (198). Thus the operator 901,) in (204) for the four cases (201) are now given by the following: _ _ _ 7, 9(u2)l(l,l)-case E 911012) = “2—39 (206) _ _ _ (1722/2 + (7 " a2) 1722) 7 9(u2)| E 622042) = 112- _ _ 20 _ 2 0 _ _ 2 _ _ _ (207) _ u; (142-7,) (u2-7,/2) - (1/2+a,) (70—17,)2+E,2/2+ (yo—172) 1722’ 90 )I e (‘) ‘ 217.42)“ (208) u 5 ll = u - , 2 (12)-case 12 2 2 8(Yo-a2)2+2(Yo-EZ) +170 0 ‘ - 0 ‘ ‘ 2,722 209 ("ZNWHW = 21%) " “2—21722+172/2+1/4Y"° ( ) In (206)-(209) the expression for 911072) is independent of 611(7) and so holds for all strain-softening materials of the type under consideration here. However the expressions for 912072), 921072) and 922072) depend on 611(7) and we have used material A given in Appendix-A for the expressions presented in (207)-(209). For 70 = 0.8 the functions 9(172): 911072), 912072) and 92,072) are shown in Figure 36. Since 70 = 0.8 < 1, from (205)4 the pure phase-II solution is not allowed and this is why 922072) is not shown. The strain 71 on first element is less than 1 whenever I 0 S 172 S N = 0.5, while the strain 72 on the second element is less than 1 whenever 70 - 1% = 0.3 S 172 S '70 = 0.8. Since the value of 172 determines the phase for each element according to Figure 35, the function 9(172) that applies is dependent on 172. For 70 = 0.8 these are: 102 (12)-case: O S 172 S 0.3, (region-B), (1,1)-case: 0.3 S 172 S 0.5, (region-A), (210) (2,1)-case: 0.5 S 172 S 0.8, (region-C). Therefore for different ranges of 172 we should follow different curves of 9(ii7). From (210) the final trace of those different functions 9072) is shown in Figure 37. Although Fig- ure 36 and Figure 37 show functions 8(1'22) for 70 = 0.8 these figures are representative for all 70 < l. The functions 9(172): 922072), 912072) and 921072) for 70 = 2 are shown in Fig- ure 38. Similarly since 70 = 2 > 1 from (205)1 the pure phase-I solution is not allowed and this is why it is not shown. The strain 71 on the first element is less than 1 whenever 0 S 172 $57 = 0.5 while the strain 72 on the second element is less than 1 whenever I 70 — N = 1.5 S 172 S70 = 2. For 70 = 2 the phase arrangements change with 172 as fol- 0.4 , .Z 7 3 I - (l,2)—case (1,1)-case (2,1)-case 992’ 72>1<+>72<171<19 71>] 0.2 — é _ 9 (t7 ) 21 2 \ 0 912(§2)\‘ -0.2 )- _ 611(H2)1 L E r “0:4 0 0.2 0.4 0.6 $2 0.8 Figure 36. The functions 9(a)) with two equal element lengths and 7,=0.8 for material A. 103 (1,2)-cas'e (1,1):case (2,1)-case 72>1‘—:>‘YZ<1 Yl] _0.4 l l l 0 0.2 0.4 0.6 (72 0.8 Figure 37. The final trace of the different functions 8(3),) for 70:08 shown in Figure 36. There is one solution to the composite trace 9(E)=0. 9L7 1 1 1 _ 1 1 1 1 1 1 m ( 2) (l,2)-casej (2,2)-case (2,1)-case 0.8 r- _ yl<1<—§>yl>l yz>l<+>yz- o -o.2 - _ - . 9213“” 922(1) _ -0.6 )- - -0.8 - _ _1 1 1 § 1 1 1 1 1 E 1 1 0 0.2 0.4 0.6 0.8 l 1.2 1.4 1.6 1.8 it. 2 2 Figure 38. The different functions 95,) with two equal element lengths and 7,=2 for material A. 104 lows: (1,2)-case: 0.. < ’2S 0. 5, (region-B), (2,2)-case: O. 5 S 172 S 1.5, (region-D), (211) (2,1)-case: 1.5 S u2 S 2, (region-C). Therefore from (21 l) for different ranges of 172 we should follow different curves of 6(172). Figure 39 shows the final composite trace of these different functions 9(172). Although Fig- ure 38 and Figure 39 show functions 8(172) for ya = 2 both figures are representative for all ya > 1. Note for a given 70 in the interval 0 S70 < 1 that there is exactly one solution, u2 = y 0/ 2, which is the pure phase solution; alternatively this pure phase-I solution can be thought of as a (12)-solution with the phase boundary at the end node s=1 or as a (2,1)- solution with the phase boundary at s=0. However for a given 70 > 1 there are three roots: one pure phase-II solution, 172 = 70/ 2, and two non-uniform solutions (see Figure 40 and- 1 007,) ' (12)-Ease; ' ' (2,2)-case ' T (2,1)-case ' o. 1- _ ylyl>l yz>l‘-;->yz<1 °"‘ ” i 7, 9 " " O N l . ,(u,)\‘/\:2/' — 1 21(“2)\‘/ 1 l -0 4 1. S = i 22(u2) S = 2- .. -0.6 e - -0.8 >- —( _1 1 1 1 1 1 1 1 1 1 0 0.2 0.4 0.6 0.8 l 1.2 1.4 1.6 1.8_ - 2 Figure 39. The final trace of the functions 9G,) for y,=2 shown in Figure 38. There are three solutions to the composite trace 9(u2)=0. 105 17 A (1.2) 7, ‘J 1 B A=(1.1) B=(1,2) C=(2,1) 1/2 (1.1) D=(2.2) A W 1/2 ‘72 Figure 40. A schematic diagram of the possible solutions 32 for a given 7, 3 l l l 7., f 7?: 0 1 7.. = rum) 2 , <\ , (2,2)-case s=1/2 is both (l,2)-case and (2,1)-case .8 ‘ ‘\ 5 = 1 (1,1)—case Yo = F103) = O 0 l l l 0 0.25 0.5: cc 0.75 0' 1 Figure 41. The possible N=2 FEM solutions as represented in (212 for 'y,=2 and for 'y,=.8 for material A. 106 Figure 41). The non-uniform solutions have the graphical interpretation of being the inter- section of the s=1/2 line in Figure 10 with the associated line of constant 70 11.3 The hard device energy function A as a function of the possible nodal displace- ments In Section 11.2 we gave a careful analysis to find the nodal displacements ii, that solve (198) for different displacements 5 and hence different 70. In this section the hard device energy function A is explored. We show directly that values of 172 that solve (198) are ex- tremum points of A, while values of that do not solve (198) are not extremum points. From (23) and the non-dimensional parameters (203), the hard device energy func- tion is given by " - _. 1 ‘72 Yo - ‘72 A042) 5 A012) = 2 (“Kl—73) + W(—1—/—2—)) - (212) Graphs of (212) for material A are given in Figure 42 and Figure 43 for the values 70 = 0.8 and ya = 2 considered in Section 11.2. Recall from (210) for ya = 0.8 that there are three possible ranges as shown in Figure 36 and these are also indicated on Figure 42. Figure 42 shows that the pure phase solution at 172 = 70/ 2 = 0.4 (Figure 37) gives the global min- imum to the hard device energy function 71. For 7, = 2, Figure 43 shows that the pure phase-II solution at 172 = 70/ 2 = 1 gives the local maximum to the hard device energy function while there are two global minima to the hard device energy functions located at the other two roots of the composite trace in Figure 39. The value 70 = l is a turning point to have either one extremum or three extrema. From Figure 37, Figure 39, Figure 42 and Figure 43 it is seen that the solutions of Section 11.2 are associated with either energy min- imum or energy maximum. The standard interpretation is that the former are stable and the latter are unstable. Thus here the pure phase-I solutions are stable (no other FEM solutions 107 1.4 — _ (12)-case (1,1)-case (2,1)-case 1.2 )- 5 .. 2 1 - ; Yl<1 Yl>l _ .e — é _ O .1: YO 3 \e. 2 = — A/ 0 0.2 - . 4 0 1 : 1 : I O 0.2 0.4 0.6 “2 0.8 Figure 42. The hard device energy functions KG?) with two equal element lengths and 7,=0.8 for material A. 1.4 L- V r ' r T‘- $ t r F r 7 q 1-2 t -* A072). - r f ‘\ - 0.8 "" 5 "' (l,2)-case (2,2)-case (2,1)-case 0.6 r- ; _. 0.4 )- _ 72>1<4+72<1 yl 71>] 0.2 *- : -4 0 0 0.12 0.14 Ofé 0.18 1 1:2 1.14 I 1.16 1.18 — 2 “2 Figure 43. The hard device energy functions KG?) with two equal element lengths and y,=2 for material A. are available), while the pure phase-II solutions are unstable. Figure 44 shows how the hard device energy functions varies with both 172 and 70. 11.4 Iteration solutions and the initial guess displacements From Section 11.3 the distributions of the hard device energy function were explored for the two special cases: 70 < l and 70 > 1. We now show how this relates to the iteration scheme given in Section 10.5. To find all of the solutions, the initial guess displacements 108 Figure 44. The hard device energy functions K62) with two equal element lengths and different 7, for material A. are chosen from the possible range 0 < 172 < yo. Consider first a displacement such that yo < 1. We then find that if the iteration scheme given in Section 10.5 is implemented, then the associated iterates “2) will follow the direction of the arrows on the energy curve in Fig- ure 42 and approach the final solution 172 = 70/ N = 70/ 2. To any initial guess displace- ment, the final solution is always the same, and is the pure phase-I solution (see Figure 42). Now consider a displacement with 70 > 1 . In this case there are three solutions (Fig- ure 39). Thus different initial guesses could converge to a different final 172. We find that the associated iterates “2)- will follow the direction of arrows on the energy curve in Figure 43. Hence there are the following conclusions for 70 > 1: (a) If the initial guess is chosen exactly at 172 = yo/N = ‘y o/ 2, .which is the pure phase- 11 solution, then it is the final solution. Otherwise for other initial guess displacements the pure phase solution cannot be found. (b) Ifthe initial guess is chosen at 172 < 70/ N = ‘y 0/ 2 , then the final solution of the (1,2)- 109 case will be found by several iteration steps. (c) Similarly if the initial guess is chosen at 172 > yo/N = 70/ 2, then the final solution of the (2,1)-case will be found by several iteration steps. In summary the domains of attraction for the initial guess to the FEM solutions are given by the well structure of the energy A. The minimum can be converged to, but the maximum requires an exact initial guess. However numerical irregularities and the finite size of e in (197) might permit an actual FEM code to stay near the maximum long enough to be found as a solution. In this section there were at most 3 FEM solutions for any 70. Thus 3 initial guesses are necessary to locate all solutions, that is all solutions of type [I I], [I II], [II I] and [II II]. However if we only need to find the possible solutions with different volume fraction, then we only have 3 general cases , and . It is found that only 2 of these can exist at the same time since (I I> and cannot both occur. Thus only 2 initial guesses would be needed to find the different volume fraction solutions. 12. Solutions of the system equation with N=3 In Section 1 1 we have explored the relations between the solutions to (198) and the possible nodal displacements for N =2. Can similar relations happen for the N =3 case? This will be explored next. Only the hard device loading is considered. Once again the energy as a func- tion of the possible nodal displacement are studied. 12.1 Formulation '5." Under a hard device loading with three elements (N =3) there are two prescribed conditions: 111 = 0, u4 = 5, and three unknown quantities to be determined: u2, u3, a, where O], E,— 1, '7le, E = E (7 ) = 671(72) (215) 2 2 2 {Ell = ’Y , 72> 1, 2 E,— 1, 13S], E3 - E303) " 011(73) >1 110 111 where “2 143—“2 5"“3 71 = T) 72 = l 1 Y3 = 1 ° (216) For each E ,- there are two choices which are determined by the strain in the i-th element. Therefore there are 2N = 23 = 8 possible phase arrangements: [1, I, I], [I], I, l], [1, II, I], [1, I, II], [I], II, I], [I], I, II], [1, H, H] and [I], II, II]. From the above observation the com- plexity of the problem is increasing very quickly as the number of elements increases. Equation (214) can be inverted and the result is given by E 8 “2 = 3 El . (217) For the pure phase case, either [1, I, I] or [I], II, II], the solution to (217) is simplified to the pure phase solution: “2 = 70|:l] = §|:1] given previously in (186). u3 2! 3 2 12.2 The hard device energy function A as a function of the possible nodal displace- ments In Section 12.1 we found that the formulation for N =3 is similar in structure to the formu- lation in Section 11.1 for N =2. However the analysis complicated by the additional degree of freedom. The possible nodal displacements that satisfy (217) will again be given by ex- trema of the energy function. Therefore we neglect a direct search for FEM solutions (as in Section 11.2 for N =2) and instead study the energy function directly. Let’s define non-di- mensional parameters in the same way as done in (203) E,- = ui/L, i = 0,1, 2, 3, (218) so that from (213) that those displacements satisfy the following relation 0l, both for material A and for material B. 1 l o 1 l I l l l l o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 17 0.8 2 Figure 47. The contour diagram of the hard device energy functions ME, 33) on the domain (219) with three equal element lengths and y,=0.8 for material A. A similar picture results for material B at ‘y,=0.8. 114 diagram. There is a pure phase-II solution to (217) at (172, 173) = (Yo/N, 27,/ N) = (2/ 3, 4/ 3) which is a hill point to the energy function. Moreover there are three valley points and three saddle points to the energy function. The maximum corresponds to phase arrangement [11, II, II]. This solution can also be thought of as a (1,2) type solution whose phase boundary is at s=0. The minirna correspond to phase arrangements [1, I, II], [1, II, I] and [II, I, I]. The first of these is a (1, 2)-type solution whose phase boundary is at s=2/3. The saddle points correspond to phase arrangements [1, II, II], [II,I, II] and [II, II, I]. The first of these is a (l, 2)-type solution whose phase boundary is at s=1/3. For any 10 > 1 the energy function and the contour diagram are similar in form to that for 70 = 2. It can be shown for material A that yo = 1 is also the turning point to have either one such critical point (minimum, maximum, saddle) or multiple critical points. In summary, for ya = 0.8 there is one solution and it is of type [I, I, I]. For 70 = 2 there are 7 solutions, and they are of type [1,1, 11], [1, II, 1], [11,1, I], [1, II, II], [II, I, II], [I], II, 1] and [I], II, H] (see Figure 49). These solutions correspond to the intersection of the horizontal lines 70 = 0.8 and 70 = 2 with the curves s=0, s=l/3, s=2/3 and s=1 of Figure 10. This is shown here in Figure 48. There is one exception to this, namely only the pure phase solution for s=l is found, the nonattainable limit for s=l at o = 6c is not represent- ed by the FEM. In all cases discussed so far (for both N =2 and N =3 cases) there was at most one solution for each phase arrangement. That is if one phase arrangement is chosen (say [1, I], 1]), then one could get the individual E ,’s in (217) by taking the appropriate formula. The resulting system equation, although still nonlinear, had at most one solution. We now show that this does not always happen. Specifically, as we now show, material B gives two solu- tions each for phase arrangements [11, I, I], [1, II, 1], [1,1, 11] when 0.97 < yo < 1 (see Figure 50). 115 3 , r 1 s = 2/ 'Y, s = 1/3 2 ,.. . H Y0 = FIKO) « . . //' Melina 1 ~ ”FEM ‘\ 1 <11.» 0 o 0.25 o.l5 = o 0.75 O' 1 Figure 48. The possible N=3 FEM solutions as represented in (2,2 for 'y,=2 and for y,=.8 for material A. Comparing this diagram to Figure 41, it is noted that now two internal constant s curves (s=1/3, s=2/3) are available for non-pure FEM solu- tions. l 1 1.6 172 2 Figure 49. The contour diagram of the hard device energy functions A632, 113) with three equal element lengths and ‘y,=2 for material A. This picture is representative for all y,>l, both for material A and for material B. 116 _ 1 1 u M A A 3 ‘(qu x m Ag,“ 8‘ D At >15 0 o A 111 0.8 - X3 0 3 A _ ° 9% .. o I - '31K* ° '.. A 0.6 - £hb,4§ ‘ a ‘5 valleypoints= 0.4 >- _. valleypoint= 0.2 ’- —1 hillpdnu= 0 l l l I o 0.2 0.4 0.6 0.8 172 Figure 50. The contour diagram of the hard device energy functions A672, 33) with three equal element lengths and y,=0.97 for material B. This is different from the contour diagram for material A with y,=0.97 (which is similar to the diagram in Figure 47). Recall from Section 5.4 for a given constant s > s that there is a minimum point trans on the curve of constant s at the point (0,,2, 7,2) . For material A we have s,,,m=2/3, while for material B we have s = 1/2. Therefore for material B the curve s=2/3 (which trans obeys s>s ) is not monotonic in 912. There exists a minimum at 7,, (s=2/3)= 2 trans s 3 (”r/3,1 -s)2/3 — + 50.9633. If 1,, <7 < 1, then there are solutions to (217) 2 2) 9:LG 2 0 with two different stresses which are not pure phase-I solutions. For example, for ya = 0.97 , there are seven solutions to (217) for material B. This is similar to what happens when ya = 2. But the main difierence is that all six multi-phase solutions to (217) corre- spond to s=2/3 in (212 and no solutions to (217) correspond to s=l/3 in (212. To show this consider yo = 0.97 for material B. The contour diagram of the energy '8\‘.. 117 s= 2/3 = = S = l 3'1“]me Ya 17(0) 0 3 valley points = = 0.25 0.5 = 0.75 O' 1 Figure 51. The possible N=3 FEM solutions to (217) as represented in (212 for 'y,=0.97 for material B. In addition to the pure phase solutions, FEM solutions are associated with the curves s=1/3 and s=2/3. There is no intersection of y,=0.97 with s=l/3 (dotted). There are 2 intersections of y,=0.97 with s=2/3, each of which corresponds to 3 phase permutations. _ "o, o o o o 0 [1511.110 x x x x x x x,,,..x-" 3.. \. “.“0 “3 o o,- o o o o o x x x x x,,...x- x o h "o o o o o x x x xxx x 1.6 [L O ' "_"‘ .11 o o o o, o o x x_ x x x x o o o o o, o X,.,.--X' "x x x x [11.1.1] 0 o o o o "x x x x x 1.2 o o o + If x x x _ o + + + + x 1"- i) + T-t- + +__ + + + «e + + + +' + + + 0-3 " JUL!) + + + + [11.1.11] - «» + + __a- + + + + + + «» + +3 + 0.4 )- -‘ cr- +’.’+ -"+ 0 [1.1.11] 1 1 l 1 0 0.4 0.8 1.2 1.6 “2 2 Figure 52. The solutions to (217) and the distributions of initial guesses in the (32, 33)- plane with three equal element lengths and y,=2 for material B. 118 0 . 8 : ¢ ¢ ¢ ¢ s t ; g :— 7; ‘r t : i t _. d» + + + + + + + + + + + + + + + + + “3 «r- + + + + + + + + + + + + + + + + .1} + + + + + + + + + + + + + + + qr + + + + + + + + + + + + + + q» + + + + + + + + + + + + + d» + + + + + o + + + + + + o + + + + + + + + + + + 4» + + + + + + + + + 0 4 <- + + + + + + + + + - w + + + + + + + + .1- + + + + + + + . + + + + + + «v + + + + + <1. I. 1) 4» + + + + u» + + + .1. + + qr + 0 J 0 0 . 4 0 . 8 ‘72 Figure 53. The solution to (217) and the distributions of initial guesses in the (32, 33)- plane with three equal element lengths and 750.8 for material B. functions MHZ, 173) is shown in Figure 50. There are also seven solutions to (217) in the diagram. The energy function at the pure phase-I solution to (217) is a valley point of en- ergy surface. There are also three other solutions to (217) that are also valley points of the energy surface. These correspond to solutions of the type <1, I, 11>. There are also three dif- ferent solutions to (217) that are hill points of energy surface. These also correspond to so- lutions of the type <1, I, 11>. Figure 51 shows these solutions as intersection points of 70 = 0.97 with the curve s=2/3 in (212. 12.3 Iterative solutions and the initial guess displacements In Section 11.4 we found that if the iteration scheme given in Section 10.5 is implemented, then the associated iterates “2. will follow the direction of arrows on the energy curves in I Figure 42 and Figure 43. In this section we implement the iteration scheme given in Section 10.5 for finding the solutions to (217) for N = 3. We find that the iteration scheme flows 119 0.8 . v . T 4 4 v . , . v v . 4 4 _ 4» + + + + . + + + + + + + + + + + + + “3 4» + + + + + + + + + + + + + + + + 4» + + + + + + + + + + + + + + + 4» + + + + + + + + + + + + + + 4» + + + + + + + + + + + + + 4» + + + + + o + + + + + + 4» + + + + + + 4» + + + + 4» + + + + + + + + + 0 4 4>- + + + + + + + + + — 4» + + + + + + + + 4» + + + + + + + 4» + + + + + + 4» + + + + 4» <1. I.I> 4» + + + + 4» + + + + + ... 0 1 0 0 . 4 _ 0 . 8 “2 Figure 54. The solution to (217) and the distributions of initial guesses in the (32, 53)- plane with three equal element lengths and 70:0.8 for material A. All initial guesses flow downhill to the <1, I, I> solution under the iteration scheme of Section 10.5. down the energy contours of Figure 47, Figure 49 and Figure 50. 12.3.1 Material A To find all of the solutions to (217), the initial guess displacements £2, {£3 are chosen around the possible ranges from (219). Under a hard device loading with 70 < 1 there is one solution to (217) only (see Figure 54). To any initial guess, the final solution to (217) is al- ways the same, it is the pure phase solution (172, 173) = (yo/N, 270/N) . We note that the iteration scheme flows down hill to the minimum, similar to the N =2 case of Figure 42. Under a hard device loading with ya 2 1 there are seven phase arrangements. Con- sider yo = 2, the (172, E3) -plane of initial guesses is shown in Figure 55. We find that we need at least 28 initial guesses (for material A) in a (172, 173) -plane to find all of the solu- tions to (217). This will be explained later. 120 7; = 0 turn 73 = 0 [11.11.11 11 72 7 = 7 c ’o. o o o o o 0 xi x x x x x x x o ‘b.’ o o o o o x x x x x.” ..x-° x o o "'o._ o o o o x x x x“ xx x 1 6 o o o "0..” o o x Xxx x x [1111] T 173 o o o o "'o_‘ 0 xxx x x x ' ' o o o o o x x x x x 1.2 o o o +___. d x x x \ .. o + + "'+[11. +1.11 + x '12 = 0 + + + + "'+ + + + e b o 8 + 411.1131] + + + + - " + + +1. + + + ‘ . 7 = 7 + + 'Iv-r + + [11.1.11] 1 3 + + + + 0.4 r- + 0+ + -‘ +. +,/ rum 3 o I 1 I 1 o ‘ 0.4 0.8 1.2 1.6 a 2 2 Figure 55. The solutions to (217) and the possible initial guesses in the (13;, 33)-plane with three equal element lengths and 70:2 for material A. The whole region is symmetric with respect to the line c—g, which is the line of 172 + 173 = 70- There is no solution to (217) in the region Aagh, r72 > 173 region, since this violates (219). There are three different regions which have their own solutions to (217). The three regions are E] aedb, D bdfc and El dehf. These regions are separated by the three lines 571, E: and if. In the whole region there are three lines which pass through these regions that have special meanings. They are 07', E: and 53. Since 70 = 2 > 1 the pure phase-1 solution to (217) is not allowed. For example, on the line 3f, the strains on the first and second ele- ments are equal, 71 = 72. Hence there are three phase arrangements to (217): [1, I, II], [11, 11, 11] and [11, 11, I] on the line 07. The distributions of the solutions to (217) on the. line H} are also similar to the two4element case. Narnely there is one initial guess which is the pure 121 phase solution, while initial guesses which are on the different sides of the pure phase so- lution will be attracted to the other solutions. Thus 3 initial guesses are needed to find all of the solutions to (217) on the line 07. On the line 571, the strains on the second and third elements are equal, 72 = 73. Hence there are three phase arrangements to (217): [1, 11, II], [11, 11, 11] and [11, I, I] on the line. It is also similar to two-element case, and three initial guesses are needed to find all of the solutions to (217) on the line b—h. On the line c_e, the strains on the first and third elements are equal, 71 = Y3- Hence there are three phase arrangements to (217): [1, II, I], [11, 11, 11] and [11,1, 11] on the line Fé which is the line of symmetry. It is also similar to two-element case which needs three ini- tial guesses to find all of the solutions on the line. In conclusion seven initial guesses are needed to find all seven solutions to (217) for a given 70 if only these three lines are con- sidered. From the above observation a three-element model which gives a two-dimensional case of initial guesses can be reduced to three two-element models which give a one—dimen- sional case of initial guesses. Point (I is a hill point and is only found from the iteration scheme if it starts from d. With the exception of this point, iterations starting on the lines 313, E: and If stay on these lines and converge to the saddle points associated with phase arrangements <1, 11, 11> while iteration starting within the open regions Cl aedb, [:1 bdfc and [:1 dehf stay on the re- gions and converge to the valley points associated with phase arrangements <1, I, 11>. But no iteration scheme will converge to a hill point and so the hill point is only found if it is the initial guess itself. There are two line segments, Zr and a}, on the line Fe. The line segment a be- longs to the region [3 bdfc while the line segment a; separates D aedb from [:1 dehf. The phase arrangement [1, II, I] is on the region [3 bdfc, while phase arrangement [11, I, 11] is on 122 the line 22’. Therefore it is easy to find the phase arrangement [1, 11, 1] to (217) since it be- longs to the region D bdfc. If the initial guess is chosen anywhere in this region then the phase arrangement [1, II, I] to (217) will be found. But it is difficult to find the phase ar- rangement [11, I, 11] to (217) since it belongs to a line. The phase arrangement [11, I, 11] to (217) cannot be found unless the initial guess falls on the line Ea. There are similar condi- tions for the other two lines, 3} and 571. Hence to find all solutions to (217) the initial guess- es are required to fall on the lines E3, 35 and a}. If we chose initial guesses in the region Aahc by forcing these initial guesses to fall on the lines (72, a) and Hf, then we found that 28 initial guesses were required to find all of the solutions to (217). It is easier to find a set of solutions which represents all the possible volume frac- tions of phase-I. This is done by neglecting the ordering the strain phases. Thus if we only consider (1,2)-type solutions, then all such solutions either involve an internal phase boundary either at S: 1/3 or 2/3, or else are pure phase solutions with phase boundaries at the ends s=0 and 3:1 (but the (1,2)-type solutions with s=0 and s=1 cannot exist for the same 70). Hence if the initial guesses are chosen on the lines a7, 1712 or E: then, at least, three different initial guesses are needed to find all of the solutions to (217) with different s. For example on the line 5}“, 71 = 72 so that the first initial guess could be (r72, r73) = (yo/N, Zyo/N) , the second one could be (yo/N —e, Zyo/N — 28) and the third one could be (yo/N + e, 270/ N + 28) where e is a small number. Note that these three initial guesses are always with E3 = 2&2 so that the strains of first and second elements are ini- tially equal. There are similar developments either for the first and third element or for the second and third element. Therefore seven initial guesses which lie on the lines a}, 372 or 22 are needed to find all solutions to (217). 123 12.3.2 Material B For a given 5 > 5mm if 70 < y¢I(s) then the analysis is similar to material A with 70 < 1, while if yo > 1 then analysis similar to material A with 70 > 1. For 7¢I(s) < 70 < 1 the whole region of (E2, E3) -initial guesses is divided into four subregions: E] aeji b, C] ehflcj, [jfcbr'k and Aijk. Since these regions: [:1 aejib, C] ehfkj, [jfcbik and Aijk give the ener- gy wells, the initial guesses which fall into one of these four regions converge to difl'erent final solutions to (217). The pure phase-I solution to (217) belongs to the triangle region Aijk. The three regions [:1 aejib, [j ehfkj and chbik are the energy wells for solutions to (217) of the type <1, I, 11>. The initial guess which lie on one of the internal boundaries (the 1 solution at point d, unless the initial guess on these internal boundaries is located exactly at one of the three hill points. These hill points are also solutions to (217) of the type <1, I, 11> (see Figure 56). o o n o o a o x f x x x x x x x x h .. D D D a cr :1 c: x x x x x x x x 3,8 .. c1 0 a c1 c: D D x x x x x x x _. .. D D D D D D x x x x x x 53 .. o 0 r: o ’13.... is .. x ’«x x x x x -- cr 0 o . cr """ x--.....x x x x 3-6 -- r: cr 0 .-. , x x ”x ""x--....::.'.',fi' — .. D o o of o 1301. xw x o o o o "of?""3;"""”‘of"""'o” vmwpdnu=¢ 1' 11> - b o o o o 9'; 9" o \x. o e 3'4 o o o o _bj' I-‘o o ‘-\ '- o o o o i? ' o \_ o o o o ,0; \"x 3.2 o o o o IIfI/ thleyPOim= _ o o o ’4' o o hillpoints=<1.1,11> a 0 0 v I l l l o 0.2 0.4 0.6 "2 0.8 Figure 56. The solutions to (217) and the distributions of initial guesses in the ('17,, 3,)- plane with three equal element lengths and 150.97 for material B. The con- tour diagram for this case is given in Figure 50. 13. FEM solutions for N>3 In this section we consider the number of different phase permutation groups described in Section 10.4. Recall from Section 11 and Section 12 that the number of phase permutation groups is dependent on yo and the number of elements N. Section 11 was devoted to the N =2 case while Section 12 was devoted to the N =3 case. In all cases, there is always ex- actly one pure phase solution. In the completely trivial one-element case, all solutions are pure phase solutions. For the two-element case, all solutions belong to one permutation group, the pure phase solution for 70 s 1, while all solutions belong to one of two phase permutation groups for 70 > 1 either <1,II> or <11,11>. For material A with a three- element case, all solutions belong to the pure phase permutation group <1,I,I> for 70 S 1 while all solutions belong to one of three phase permutation groups for 70> 1, either <1,I,11>, or <11,11,11>. We also recall that the solution for N =3 with material B was more complicated. In this section we will explore the N >3 cases to see whether there are similar results. The number of elements determines the number of nodes and hence the possible phase boundary locations. Thus the number of elements gives the available constant s curves in 912 according to (193). There are N + 1 phase permutation groups (see Figure 34). Each of which has C(N, m) members. The intersections of constant 70 lines with con- stant 3 curves in $212 determine the FEM solutions. Each such intersection with a constant 5 curve corresponds to exactly one solution of type (1, 2). However it also corresponds to the other C(N, m)-1 members of the permutation group which are not type (1,2) solutions (where s = m/N c: m = 3N). If N is large enough, then some of the curves of constant 5 curves will obey s > s and certain lines of constant 70 will intersect the same curve of constant s twice. Iran: 124 125 For example this will occur for 3: 1-1/N if 1-1/N>s > 1-2/N. Thus the condition for trans such a double intersection to occur is that N obeys N> —l—. (221) l — strans For N=3, this occurred for material B (s = 1/ 2), but did not occur for material A trans = 2/ 3). This also shows that the inequality in (221) is strict. (Strans For material A (s =2/3), condition (221) is first obeyed when N =4. Then the trans curve of constant 3: 3/4 satisfies the condition 3 > s and some range of 70 will inter- trans" sect this curve twice. Note that each such intersection represents C (4, 3) =6 solutions, and both intersections give solutions from the same permutation group: <11111>. We shall say that these two intersections give two solutions that are distinct to within a permutation (per- mutationally distinct). From Appendix-A we have Y¢I(3 / 4) < 1. Hence we conclude that the number of solutions that are distinct to within a permutation is: exactly one for 10 < 7¢I(3 /4); two for ya = 7¢I(3 / 4); three for 7¢I(3/ 4) < 70 < 1; two for 70:1; four for ya > 1 . It is easy to show the similar development for a higher number of elements. The per- mutationally distinct solutions for N from N =1 to N =8 with either material A and with ma- terial B under a hard device loading are shown in Figure 57. The number inside the L] in- dicates the number of the permutationally distinct solutions. Note that the positions of 7o 2 are not in scale. The number of permutationally distinct solutions in Figure 57 is found by counting the number of intersections that a constant 70 line in (212 has with the internal constant 5- curves: 1/N, 2/N,..., (N -1)/N that are interior to (212. This number is then increased by one to account for the single pure phase solution. Consider, for example, the eight-element case with material A. There are two values of s, 7/ 8 and 6/ 8 , which are greater than the value of s = 2/ 3. Hence these two constant s-curves have one intersection point, which is trans 126 =8 N: N=6 1v: =2 To 3 0.5 To I l . I" m - £331 material A . lll ‘ 70,0") [3) ‘7”(6/8) [5) ‘ [a] , ' 7. material B Ill 170,0") [31 1.1(6/8) [SI 1715‘” 3) [71 17.1““) (a) [21 [4) l6] ‘91 . . ‘2' m [:31 material A z m ‘ 75“”) 13) I 74’”) 151 m '———*—__9 7. material B [11 17.5“" I31 1 701”” 151 m'or‘m’ "1 . m . I [21 :3] mammal A l" ‘ 7.55/5) {3] 1 1.14/6) [5] 7. . f 4/ material B [11 [21 7.115") [3) "'11.; 6) [5) IIII 1.53/5) [6] . [21 I, (21 mammal A - “l ‘1‘: S) [3] [5] *.—p 7. material B - m 1’4”” I31 14E "(W ‘5’ : m [2) [31 , [21 [521 material A . Ill £1.10“) [3) . [41 _ ¥ 7. material B m [inf/0 [31 1.10") (4] 121 . in material A “1 1.52/3) [3) 1. material B m 1 1.5”” I31 . m A 121 1,21 . m material A ‘ m J, 121 7 . f 1 (v2) a . mammal B “I , '1 [2) m material A m m 4 '- material B I” m Figure 57. The permutationally distinct solutions for N from N=1 to N=8 with material A and material B under the hard device loading. The number inside the L] indi- cates the number of permutationally distinct solutions. Note that the positions of You are not in scale. The T give values of 70 at which the number of permu- tationally distinct solutions suddenly changes. Often the number of permuta- tionally distinct solutions exactly at the T values is different from the number of permutationally distinct solutions that is available on either side of the T values. 127 not the end point, with the horizontal line 70 = 1. The number of the permutationally dis- tinct solutions increases from 1 to 8 as 70 increases from 0 to values greater than 1. For 70 < y¢2(7/ 8) there is no intersection point with internal constants curves and there is only the pure phase-I solution. When 70 = y¢I(7/ 8) , the curve of constant s= 7/8 is just tangent to the horizontal line of constant 70 and there are two permutationally distinct solutions, the intersection point and the pure phase-I solution. Since 7¢I(7/ 8) < 7¢I(6/8), if 70 is in the range of 7¢I(7/ 8) < 70 < y¢I(6/8), then the curve of constant s= 7/8 intersects the hori- zontal line of constant 70 twice and there are three permutationally distinct solutions, two intersection points with s=7/8 and one pure phase-I solution. When the value of 70 reaches to 70 = 74,2(6/8), the curve of constant 3: 6/8 is just tangent to the horizontal line of con- stant yo and there are four permutationally distinct solutions which are given by the two previous intersection points (s=7/8), one tangent point (s=6/8) and one pure phase-1 solu- tion. Similarly when yd, (6/ 8) < 70 < 1 there are five permutationally distinct solutions 2 which are given by the four internal intersection points (two for s=7/8 and two for s=6/8) and the one pure phase-1 solution. The load value 70 = l is more interesting, since then half of the internal intersection points now coincide with the pure phase solution. Thus it is difficult to distinguish whether the solution to (217) is a pure phase-I solution, a pure phase- 11 solution or a mixed state solution. Hence there are three permutationally distinct solu- tions which are given by the two remaining internal intersection points (s=7/8 and s=6/8) and the one “pure phase-I” solution. When 70 > 1 it is similar to a general number of ele- ments case and there are eight permutationally distinct solutions. These eight permutation- ally distinct solutions are given by the seven internal intersection points of the seven differ- ent curves of constant 5 with the horizontal line of constant 70, along with the one pure 128 phase-11 solution. The discussion above shows the maximum possible number of permuta- tionally distinct solutions. From the FEM point of view, the possible permutation groups that are presented for each N depend crucially on yo (see Figure 57). 14. Numerical results using the finite element method for material A In this section several numerical examples using the finite element method are studied. The goal is to use the two separate criteria, the minimum energy solution criterion and the ki- netic relation, to pick out an unique solution as done in Sections 6, 7, 8 and 9 for the con- tinuous problem. For both criteria and for a given loading all of the permutationally distinct solutions which are explained in Section 13 are considered as found and so are available as candidates for satisfying the criteria under consideration. The selected solution is chosen from all of these permutationally distinct solutions so as to best satisfy the two separate cri- teria under consideration. The main issue in this section is how to implement these two cri- teria numerically. Recall from Section 13 that we only consider the permutationally distinct solution as given in Section 10.4 and so neglect the ordering of the phases. Thus we only consider pure solutions and (1,2) type solutions so that at most only a single phase boundary is present. In this section a quasi-static motion is considered and so we introduce discrete time steps: t1, t2, . .. For such quasi-static motions the phase boundary is allowed to change its location. Let a discrete time step begins at t = ’1' and ends at t” 1. At both t,- and ti“ the values of s are restricted to the set of nodal values 1 3 N}, (222) for a N -element model with equal element lengths. From the given conditions, either initial conditions or a calculation from the previous time step, the FEM solution at time t,- is re- garded as known. Thus s(ti) e m is also known. It is desired to determine the FEM solution at t = ‘1' + 1. By using FEM we seek to find s(ti +1) 6 m which satisfies the criterion under consideration in the best possible way. 129 130 14.] Minimum energy criterion For the minimum energy criterion the solution path is independent of the past history of the system. The associated values of the energy A are computed for all of the possible phase arrangements __ 1 N un-t-l—un Ma], 14,, uN+1) _ N Err—ET). (223) n: The minimum energy criterion gives the solution which has the least energy from among all of the available solutions at each instant time t. For example let’s consider a hard device loading with a monotonically increasing end elongation. For the two-element case the min- imum energy path first follows the pure phase-I solution, the line of s=1 in (212, from the origin 0, (o, 70) = (O, O), to the point P, (0, yo) =(1,l), secondly it traces the curve S: 1/ 2 in the (6, ya) -plane (Figure 58). In particular there is no jump in the solution path. This is 2 r— 1 1 r ‘jir ' r l r Ya =0 F \ A’ 3:0 s=1/8 . II - =8 7 1.2 _ s=6/8 ii 0 N=20 3—2/32 \\ s=3/8- 1 '- S=7/8 } 7 g - ~__~ .._-i' A1 0 8 " i :E .B 4 F o 6 L s_1.__}'§ C: ‘ ‘.— 3:1 0 4 - g + 0-2 " -1 0 0 r r r r i. r r r r 0 01 02 0.3 04 0.5 06 O7 08 09 1 Figure 58. The solution paths in 9,2 obeying the minimum energy solution criterion for the two-element, the four-element, the eight-element and the twenty-element cases under a monotonically increasing hard device loading with material A. 131 because the available FEM solutions are only the pure phase solution and multi-phase so- lution on the curve s=1/2. For the four-element case the available FEM solutions are the pure phase solution as well as the curves s=1/4, s=l/2, s=3/4. Recall from Section 6 that for 70 > 1 the energy A decreases as 0 decreases on the lines of constant 70. Thus for ya > 1 the lowest energy path will be on the leftmost available s-curve (the path closest to o = cc). However for cc < yo < l the energy is minimum as o -> cc (the unattainable limit), with the energy de- creasing toward this limit only for CC 1, it is not always the latter (namely s = l — l/ N ). Instead a transition will take place at ya = 7¢I( l - l/N). Before this transition the minimum en- ergy is given by the pure phase solution, while after this transition the minimum energy path will be the leftmost available path (the path s=1-1/N). For N =4 the path first follows the curve 3:] from the origin as in the two-element case. However at the point A the value of A is the same as the value of A at the point D (which is the intersection point of the curve Q3 and the curve s=1 — l/N = 3/4). Thus the path will jump from the point A to the point D and from then on will follow the curve s=3/ 4 as 70 continues to increase. For the eight-element case the minimum energy path first follows the curve s=1 from the origin as in the two-element case and the four-element case. At the point B the value of A is the same as the value of A at the point E (which is the intersection point of the curve $3 and the curve given by S: 7/ 8). At this value of 70 the path will jump from point B to point B and finally follow the curve s=7/ 8 as 70 continues to increase. ' Accordingly if the number of elements is increased, then the path will approach the 132 unattainable limit 6c: 1 / 2 which matches the conclusion of section 6.5 as we expect. For example for a twenty-element case, the minimum energy path will be even further to the left. The path first follows the curve s=1 from the origin as in the previous cases. At the point C the value of A is the same the value of A at the point F (which is the intersection point of curve E3 and the curve given by s=19/20). At this value 70 the path will jump from the point C to the point F and finally will follow the curve given by s=19/20 as 70 continues to increase. 14.2 Kinetic criterion For the kinetic criterion the solution path is dependent on the rate of loading and so the FEM implementation requires more work. Recall from Section 12 that the kinetic criterion gives (129). Thus if the kinetic criterion is used as the extra condition to find the unique solution from the available FEM solutions, then this solution should be chosen so as to best satisfy (129) from the available FEM solutions. Recall from Figure 25 for the (1,1)-case, (1,2)-case and (2,2)-case combination that the value of 3 always increases and this increase is continuous, except for possible sudden changes at o = 1.For example if the initial value 3 = s(t,-) = 1 and 0' at 1 from t,- to t” 1, then the only possible solution at the next time step is the solution with the final 3 = s(ti+1) = l.Butifthe initialvalues = s(t‘-) < 1 and 0:1 from ti to ti+1,thenfrom (140) the solution at the next time step should have the phase boundary 3 = s(ti+1) at the location SUI-+1) = s(t,.) + k j j: *1 f12(0(t))dt. (224) Since the FEM is a discrete representation of a continuous physical system, the FEM pro- vides the several discrete values of s which are given by (222). Equation (224) can be ap- proximated by many possible ways in a discrete fashion. Given an FEM solution at t = ti, 133 so that 700,-) , o(ti) and s(ti) are known, and given 70(‘1' + 1), the problem is to determine the FEM solution at t = t,- +1 whose values GUI-+1) and s(t‘- + 1) best satisfy (224). 14.2.1 Several algorithm to find s (rm) One way to implement (224) is by the following algorithm: Algorithm 1a: (a) Compute s by est s... = so.) + lemma.» (1,, 1 — n) . (b) For each available FEM solution at t=ti + 1 extract s(t‘- + 1). (c) For each result in (b) compute “1 = 1501+ 1) " sestl' (d) After cycling through (b) and (c), select the FEM solution that minimizes ‘1’. Physically, this algorithm estimates the integral (224) by assuming that the stress is con- stant from t,- to t Hence the phase boundary moves with constant velocity which is the H 1' initial velocity at t= ‘1" Note that Algorithm 1a is equivalent to the following: choose the FEM solution at t= t,- + 1 that minimizes ‘1’”, where ‘1’” is when... 1), so), on). r, . 1 - n) = |s(r,. 1) - s0.) - k 012m,» (t,, 1 - t.) >|. (225) In Algorithm 1a the stress is assumed to be the constant which is the stress at t= ti. How- ever the stresses of the available FEM candidates at time t=ti +1 may be different. There- fore Algorithm ]a can be modified by assuming that the stress is equal to GUI-+1) in the integral (224). This gives the following algorithm. Algorithm 1b: (a) For each available FEM solution at t: ti + 1 extract s(t,- + 1) and 0(ti + 1). 134 (b) For each choice in (a) compute s by est Sest = S no for material A and s' —) 5k/ 8 for material B when 0’ -) 1/ 2. Therefore the phase boundary will start to move slowly but will finally move very fast. 139 3 .4 _. . , r r l r .1 3 .2 - \ 0 ...... exact path 4 3 N: .. 7" r \ A....Algorithm 1a 2 .8 *- .. 2.6 L « X.....Algorithm 1b .. 2 " " *.....Algorithm 2 ‘ 2 .2 '- _. 2 ,. CI...Algorithm3 g 1 .8 '- _ 1 .6 '- 1 1 .4 - _ 1 .2 r- _. 1 I 1 ~ 0.5 0.6 0.7 0.8 0.9 0’ 1 Figure 60. The solution paths obeying the kinetic criterion for difierent algorithms under a hard device loading with 10 loading time steps. Here the conditions are given by (229) and (231) with material A and k=4.05. To see the effect of the number of elements on the FEM process, three different number of elements: N = 2, 4 and 8 are considered. Then for each different N, the influence of changing the number of time steps is examined by considering 10, 20 and 30 time steps (see Figure 61). In Figure 61 we have used a Combination Algorithm which is a combination of Al- gorithm 2 and Algorithm 3. In the Combination Algorithm at the odd numbered time steps the integral in (224) is computed using Algorithm 2 while at the even numbered time steps the integral in (224) is computed using Algoriflrm 3. For N =2, the phase boundary location only has three choices, 0, 1/2 and 1. Figure 61(a) shows that s increases monotonically. At the beginning the solution path will follow the curve s=0, jump to the curve s=1/2, and then follow the curve s=1/‘2. Due to the avail- ability of only one internal s curve for a two-element discretization there are big errors be- tween the FEM solution paths and the exact solution path. Therefore increasing the number 140 of time steps gives little improvement for N =2. But for N =4, there are five choices for the phase boundary location s: O, 1/4, 2/4, 3/4, and 1. Moreover s also increases monotonically. The solution path will first follow the curve s=0, jump to the curve s=l/4, and then follow the curve s=1/4. After that the solution path will jump to the curve s=2/4 and so on. The errors between the FEM solution paths and the exact solution path have been reduced. Obviously increasing the number of time steps gives more improvement for N =4 than it does for N =2 (Figure 61(b)). For N =8, there are nine choices for s: O, 1/8, 2/8,..., 7/8, and 1. Again s increases monotonically. The solution path will first follow the curve s=0, jump to the curve s=1/8, and then follow the curve s=1/8. After that the solution path will follow the curve s=2/8 and so on. The solution paths for N =8 are similar to those for N =4 (Figure 61(c)). Figure 61 shows that improvements take place with an increasing number of time steps when the number of elements is held fixed. Three time step discretizations for N =2 are shown in (a); those for N =4 are shown in (b); those for N =8 are shown in (c). Corre- spondingly if the nine FEM solution paths in Figure 61 are rearranged, then the improve- ments obtained by increasing the number of elements N when the number of time steps is held fixed are shown in Figure 62. 14.2.4 Three ramp loadings with different loading rates In the previous subsection, the finite element program is checked against an example which has an “exact” solution. The results are very close. The influences of the loading rates on the solution paths as computed by FEM are explored in this section. Let’s consider three different ramp loadings with different loading rates a which are similar to those loading as given in Section 9.3 and shown in Figure 63. Our goal is to show that FEM gives the rate dependence predicted in Section 9.3. Again the Combination Algorithm which alternates between Algorithm 2 and Al- 141 T 0 ......... exact path A 10 timesteps 20 timesteps l 1 l l l I I I 1 1 I I I I I 1 1 HNb®mNNfi0~mwa o m o a O \t A N v 0 (D O ITTIIITTIIT O U‘ C 0‘ O J 3 0 co O 0 Q ..o NNNMwa \ ITIIIIIITII HHHH HNha‘mNNbO‘mh’Nfi O .5 0.6 0.7 (C) 0.8 0.9 o» 1 Figure 61. The efiect of refining the spatial discretization. The solution paths obeying the kinetic criterion for the two-element (a), the four-element (b) and the eight-ele- ment cases(c). This is done by a hard device loading with three difl’erent load- ing time steps. Here the conditions are given by (229) and (231) with material A and k=4.05. The curves indicated by are available for the value of N being considered. The combination algorithm was used for implementing the kinetic criterion. 142 3.4 » ' ‘ .. 3.2 ~ _ 0 ....exact path 4 yo 3 _ lOtrmesteps g 2_3 _ + ”"N=2 _ 2.6 >- .. 2.4 _ A ..... N=4 .. 2.2 - _ 2 h- =8 _ 1.8 e - 1.6 ~ ... 1.4 ~ .- 1.2 r- 10.5 0 6 0 7 (a) 0 8 0_ 0’ 1 3.4 T 3.2 2 ‘Yo 3 - 2.8 _. 2.6 .. 2.4 1 2.2 a 2 _ 1.8 _ 1.6 _ 1.4 _ 1.2 - 1 O. O- 1 3.4 _ 3.2 - To 3 ~ 2.3 4 2.6 .- 2.4 _ 2.2 _ 2 d 1.3 4 1.6 3 1.4 _ 1.2 - 1 __ o. a 1 Figure 62. The effect of refining the time discretization. The solution paths obeying kinetic criterion for the two-element, the four-element and the eight-element cases. This is done by a hard device loading with three different loading time steps, 10 time steps (a), 20 time steps (b) and 30 time steps (c) for initial con- ditions given by (229) with material A and k=4.05. The combination algo- rithm was used for implementing the kinetic criterion. 143 Figure 63. Three loadings with different loading rates. gorithm 3 is used in the finite element program. The following conditions are assumed: 70(0) = 0, 0(0) = 0, 5:0 when 0:1. (235) We consider material A and k=lO. The values of a are given by 2,1 and 2/3 in loading 1, 2 and 3 respectively. Figure 64 shows the FEM solution paths obeying the kinetic criterion (129) for a two-element discretization (a), a four-element discretization (b) and an eight- element discretization (c) under these loadings. This figure shows that the solution path of the fast loading rate, loading 1, will follow the curve of pure phase solution s=0 longer than those of the lower rates and always remains to the right of the solution paths associated with lower loading rates. Similarly the solution path of the slowest loading rate, loading 3, will follow the curve of pure phase solution s=0 shorter than those of the high rates and always remains to the left of the solution paths associated with high loading rates. Note also that the loading 2 (which is located between loading 1 and loading 3 in Figure 63) has a solution path in Figure 64 that is also located between those of loading 1 and loading 3. Figure 65 shows the effect of the number of elements under the same loading. The solution path associated with a larger number of elements always leaves the pure phase so- lution curve s=0 earlier and has more variation. This is because these solution paths have ‘1 O-‘HHH oooo ouaompmammu 5‘ HHD—‘H oubamewt-ocom OOOO rrrrs‘ oooo owhmmHNAo‘mN + ....loading 1 A loading 2 x ....loading 3 0 .12 1 .1 o .‘2 t ‘4 (b) 0 .‘6 0:8 0' 1 . _ r r ,r r r m , N=8 : h 0 .‘2 0 .14 (C) 0 .‘6 0:8 0' 1 Figure 64. The efiect of refining the spatial discretization. The solution paths obeying the kinetic criterion for the two-element (a), the four-element (b) and the eight- element cases (c) under the three hard device ramp loadings. Each loading has a different loading rate for conditions given by (235) with material A and k=10. The curves indicated by J. are available for the different element discret- izations. The combination algorithm was used for implementing the kinetic criterion. 145 loading 1 + N=2 A ..... N=4 0.5 0.6 0.7 (C) 0.8 0.9 0' 1 Figure 65. The solution paths previously displayed in Figure 64 are here displayed (for “(>1 and material A). Here we group the FEM paths associated with the same loading rate on the same frame. Within each frame we see the effect of increasing the number of elements. 146 more choices to adjust themselves to match the requirement of approximating (224). As a consequence it will give a better result. The exact solution would be expected to leave the curve s=0 immediately. The solution path for the lowest rate will be closer to the curve s=1, which is the unattainable limit. This means that the solution paths will also approach the minimum en- ergy solution while the loading rate is very low. But there is a difference, namely the solu- tion path obeying the kinetic criterion cannot jump away from a pure phase solution until it passes the yield point (6, yo) = ( 1, 1) while the solution path obeying the minimum en- ergy could jump earlier. 15. Numerical resuns using the finite element method for material B In the previous section several numerical examples for material A were studied. Several similar numerical examples for material B are presented next. We find that the results are similar. 15.] Minimum energy criterion Let’s consider the loading as given in Section 14.1 under the minimum energy criterion. We find that the solution paths for material B, although similar to those of material A, will involve an earlier jump (see Figure 66)). 15.2 Kinetic criterion For the kinetic criterion, the main difference between material A and material B is that functions for the driving tractions are different. For material B we have s = 0 when 0 = 1 2 . Yo I I I I II' I r I 1.8 - + N= 3:4”; I; :0 - '57 s— /8 1'6 '- A ..... N=4 s=5/8 ‘ 3:23 - ' =8 . S= 1.2 - s .3 . _J 0 N=20 =6”; 4?» - \ P s=7/8 ' - " A o a _ ._ M 4 :- B o o »- ‘121'4‘ - s'l—V' C ‘— s=1 0.4 4 _ E 0 2 "O -‘ 0 l J l I 1' l I l l o 01 02 0.3 04 0.5 06 07 08 ()9 1 Figure 66. The solution paths in 0,2 obeying the minimum energy solution criterion for the two-element, the four-element, the eight-element and the twenty-element cases under a monotonically increasing hard device loading with material B. 147 148 and we have s' —-) 5k/ 8 when 0 -) l/ 2. Note that, therefore, the velocity of phase bound- ary for material B is always finite while the velocity of phase boundary for material A could approach infinity. Hence FEM solution paths for material B will typically stay longer on the curve s=0 of pure phase-11 solutions. 15.2.1 Computation of an “exact” solution path (o(t), 70(0) which satisfies the kinetic criterion Since material B is different from material A, the same applied stress 0(t) will give a dif- ferent kinetic response 70(t). Thus if o(t) is again given by (231) and Figure 59(a), the function 70(t) will no longer be given by Figure 59(b). Recall from (143) that we find that 70(t) is given by 700) = (0(1)2 + 1) (236) r,,(o(t)) + (1‘,(o(t)) - 1‘”(o(t))) (4.05 I (1, {—2—— - ./(20(t) - 1) )dt). Recall from (b. 1), (b4) and (b.24) for material B that the driving traction Rd) is given by 2 from) = 5102—111 - mom - 1). (237) Therefore from (140) and (237) the location of phase boundary is given by 1(1) = 5(0) + kjgiroondt = 4.05j3i1o(r))dt 2 = 4.05]; (Ex—l) — ((20(t)-1))dt. (238) To find the integral in (236) Simpson’s rule is again chosen with three different number of time steps,10, 20 and 30, as t increases from 0 to 1. As in Section 14.2.2, 0(t) is given by (231) and k=4.05. Figure 67 shows these three results for 70(t) , they are all very close. But the corresponding response 70(t) of material B is smaller than that for material A (Figure 59(b)). 149 0 o 0T2 of; (a) 0:6 038 t 70(1) 111111 HHHH l l llllllllll 0000 1111 owaamewaameN o 0.1 0.2 0.3 032 of: 036 0.‘7 1138 o.‘9 (b) t Figure 67. (a) The loading o(t) given by (231) which is also considered for material A. (b) The corresponding response function 70(t) which is given by (236) using (231) for material B with k=4.05. Three different total number of time steps, 10, 20 and 30, are considered in the numerical computation of this “exact” solution. 15.2.2 Solutions obeying the kinetic relation obtained by using the finite element method Let’s consider three different number of elements, N =2, 4, 8, and three different number of time steps, 10, 20 and 30, as given in Section 14 for material A. In what follows the Com- bination Algorithm is also considered. The FEM solution paths of these cases for material B are similar to the corresponding FEM solution paths for material A as shown in Figure 68. Figure 68 shows the effect of refining the spatial discretization. The phase boundary lo- cation of the exact path for material B is distributed form s=0 to s=3/8 (while those for ma- terial A ranged from s=0 to s=7/8 (Figure 60)). Correspondingly Figure 69 shows the ef- fect of refining the time discretization. 150 2.2 )- \\ -( ‘ \\ S=0 =2 0 ......... exact path ‘ .. A 10 time steps x 20 time steps + 30 timesteps 0.5 0.6 0.7 (c) 0:8 0.9 ’- 0- 1 Figure 68. The effect of refining the spatial discretization. The solution paths obeying the kinetic criterion for the two-element (a), the four-element (b) and the eight- element cases(c). This is done by a hard device loading with three difierent loading time steps. Here the conditions are given by (229) and (231) with material B and k=4.05. The curves indicated by '1' are available for the value of N being considered. The combination algorithm was used for implementing the kinetic criterion. 151 2 -2 0 ...exact path - Va 2 + N=2 ‘ 1.8 A ..... N=4 .4 x N=8 H 'a. 1 ... N 1 Figure 69. The effect of refining the time discretization. The solution paths obeying the kinetic criterion for the two-element, the four-element and the eight-element cases. This is done by a hard device loading with three different loading time steps, 10 time steps (a), 20 time steps (b) and 30 time steps (c) for initial con- ditions given by (229) with material B and k=4.05. The combination algo- rithm was used for implementing the kinetic criterion. 152 15.2.3 Three ramp loadings with different loading rates We saw in Section 14.2.4 that the loading rates do have different effects on the solution paths of material A. Naturally it is desired to see whether the loading rates have a similar effect on material B. Let’s again consider the three different loadings with different loading rates as given in Figure 63. Figure 70 shows the effect of refining the spatial discretization. Figure 70 shows the solution paths obeying the kinetic criterion for the two-element (a), the four-element (b) and the eight-element cases (c) under the three hard device ramp loadings. Each loading has a different loading rate for conditions given by (235) with material B and k=10. Figure 71 groups the FEM paths associated with the same loading rate on the same figure. Within each frame we see the effect of increasing the number of elements. The solution paths for mate- rial B are also similar to those for material A. 153 70le + ...loadingl = g : 1-6 ’ A....loading2 ‘ 1.4 4 32‘ 4 1.2 ,_ x ....loading3 f 3 1 4 0.8 4 i - 0.6 4 f s=7/8 4 0.4 4 0.2 4 _ 0 0 0.11 0.12 0.13 0. 4 (a) 01 g 04.7 0 8 0.190. 1 1101 To 2 " " 1.8 4 4 1.. _ N=4 4 1.4 4 _ 1.2 4 .1 1 ~ I 0.8 4 _ 0.6 4 4 0.4 4 4 0.2 4 4 0 0 0.11 0.12 0.13 0. -4(b) 0:. 5 016 0.17 (is 0.19 O' 1 “Hill 7, 2 ‘ 1.8 4 _ 1.6 4 N: 4 1.4 4 .. 1.2 4 . 4 1 4 . 0.8 4 5 4 0.6 4 4 0.4 4 4 0.2 4 Q 4 00 011 0.2 0.3 0.14 (c) 015 0.16 0.17 0:8 0:9 0 1 Figure 70. The effect of refining the spatial discretization. The solution paths obeying the kinetic criterion for the two-element (a), the four-element (b) and the eight- element cases (c) under the three hard device ramp loadings. Each loading has a different loading rate for conditions given by (235) with material _B and k=10. The curves indicated by 1' are available for the different element discret- ization. The combination algorithm was used for implementing the kinetic cri- terion. 154 0.5 0.6 0.17 (b) 0.8 0.9 a 1 0.5 0.6 0.7 (C) 0.8 0.9 O. 1 Figure 71. The solution paths previously displayed in Figure 70 are here displayed (for y>l and material B). Here we group the FEM paths associated with the same loading rate on the same frame. “Within each frame we see the effect of increasing the number of elements. 16. Conclusions and recommendations for future work 16.1 Conclusions From this study we may draw several conclusions: 1. The nonuniqueness of the analytical solutions to the boundary value problem trans- lates into nonuniqueness of the numerical FEM solutions to the boundary value prob- lem in the simplest model (a bar). 2. It is possible to successfully implement two separate selection criteria, either the ki- netic relation criterion or the minimum energy solution criterion, to the phase transi- tion problem by using FEM. 3. If the number of elements is increased, then the number of internal s curves in (212 also increased. Therefore the FEM solution paths are more smooth. Finally the solu- tion path will approach to the exact path as N gets larger. But at the same time the number of initial guess solutions needed to obtain the available FEM solutions in- creases dramatically as N gets larger. 4. Even though in this study two separate criteria are implemented on the FEM process to a one-dimensional problem, it is reasonable to apply these two separate criteria to other phase transformation problems using FEM. 5. In this study we only considered quasi-static motions which neglect the effect of in- ertia. As shown in Lin and Pence [1993LL] these should give the large time solution of the dynamical problem. Thus quasi-static FEM might give good predictions for the ultimate fate of fully dynamical solutions. 16.2 Recommendations for future work A lot of issues can be explored after this study. For example these could include the follow- ing: . 1. In this study the hard device problem is the main concern. However it would be inter- 155 156 esting to explore further the soft device problem. 2. In this study a bar which has unit cross-section area is explored. Therefore the stresses on each element of the bar are equal in the FEM analysis. But if a tapered bar is con- sidered by the similar process described in this study then one should have N different stresses in the N -element case. Therefore there could be N different solution paths which are obtained from each elements. 3. In this study a linear relation between the phase boundary velocity and the driving traction is considered. There are other kinetic relations to explore the strain-softening behavior of the bar. Moreover one can try other criteria. 4. The boundary value problem given in Figure 3 are solved by quasi-static motions. However it might be possible to solve this kind of problem by considering fully dy- namical motions. 5. It might also be possible to find a spring-dashpot model to represent the strain-soft- ening behavior for both the hard device and the soft device. 6. In this study two example materials which can resist an infinite strain are concerned. But it is possible to consider hypothetical example materials with the specific material behavior which model different processes. For example choosing a material which only can resist a finite strain could model damage accumulation within a material. Appendices Appendix-A: Material A In this Appendix and the next appendix, two example materials which have the behavior of strain-softenin g are introduced. This simply means that particular functions for the function 611(7) of (5) are considered. Material A is taken to be the following: 5,0!) = Yr 751, 5(7) = _ 1 1 (a.1) 611(7) = 2+7? y>l. Note that the constant CC = 1/ 2 so that the behavior of strain-softening 611(7) —> 1/ 2 when 7 —) 0°. The first derivative of 5” (y) is given by 6,,’ (r) = i < o, (42) 212 which shows the strain-softening. The second derivative of 611(1) is given by 6””(7) = $>0, (a.3) which shows the curvature restriction. From (a. 1)2 the inverse function of 611(7) is given by F o) - ——1— (a 4) ”( — 20 - l ' ° From (7) and (a. l)2 the energy density function W (7) for y) 1 is given by 1M0 W(y) = 2+T’ ify>l. (a.5) Thus from (a.4) and (a.5) the energy density function W (1‘1,(0)) can be rewritten as l l l . 5 (m +1"(-2?:—1)), While r"(0') > 1. (3.6) W (l‘,,(o)) = Note from (10), (a.4) and (a.5) that the limit We is wee lim (Wm —yo() = Iim é—lnm = 4». (a.7) re“ ra~ 157 158 This means that the area under the strain-stress curve that is also above the horizontal line 0 = 0c is infinite. According to (a.4), the equation (23) which involves the boundary con- ditions u(O) = 0 and u(L) = 50 can be rewritten as 05+0(1-s) = 0 = yo, (1,1)-case, l 0s+ 201—1 (1 —s) — yo, (1,2)-case, —-——1 s+0(1 -s) - (21)-case (a.8) '7'? 20-1 - 70’ r 9 ——1 s+—l—(l-s) --1—- (22)-case 20-1 20-1 - 20-1— 70’ ’ ° t From (a.8) it is found that (20- 1)n-l .. = iv 312‘0’70’ (20+1)(o-1)’ ( 9) . (20-1)0-—yo(20-l) 3‘ s21(0, 7.) = (20+ 1) (0 — l) ’ while 311(0, 70) and 322(0, 70) do not have meaning. Alternatively (a.9) can be found di- rectly from (25) and (a.4) using F,(0) = 0. From (24), (a.4) and (a.6) the hard device energy function A for material A is given by o2 . All = —2' : A11(09 YO)’ (1,1)’CaSC, l 1 2 -———+ln(———) 0 20-1 20-1 .. A12 = 75+ [ 2 ](1—s) = A12(0,'yo), (1,2)-case, 1 +1 ( 1 ) (a.10) _ n _— 20-1 20-1 02 .. A21 = [ 2 ]s+-2—(1—s) = A21(0,'yo), (2,1)-case, l 1 ————+ln(—) 22 = 2 = A22(0r70)r (2,2)-case, where s is given by (a.9). The soft device energy function E for material A can be obtained 159 from (26) and (a. 10) 511 = A11(0,yo)-0yo, (1,1)-case, 312 = A12(0,yo) - 070, (1,2)-case, .. (all) 321 = A21(Gr Yo) - 0709 (2’1)'Cascr 322 = A22(0, yo) - 070, (2,2)-case. Note that the difference between the hard device energy function A and the soft device en- ergy function E is 070. A.l The solution region (212 for material A After material properties and the energy functions being introduced, from (29) and (a.1) the solution region 012 is given by 1 . 1 1 1 . £212: {(0,yo)|§<0 s trans trans the location of the minimum point on the curve of constant s can be found by 2(1—s) S— __2 = 0, (a.13) (20— l) where the root of (a.13) is given by _ _ 1 2 (1 -s) 0 — 0¢I(s) — 2 (1+ s ). (a.14) Hence from (47), (a.4) and (a.14) the minimum 10 on the constant s curve is given by Vol, = §+J2s(1-s) E114,20). - (a.15) Therefore from (a. 14) and (a. 15) the minimum point on the curve of constant s > sum, rs given by 160 1 1 2 1-s (0.1,) = (0,50,11,50) = (2421—13—3) 3+ +../2s(1—s)) (a.16) Note from (a.14) that 0 = 0102 = 1 when s = 3mm = 2/ 3. Note also that the stress in (212 must obey 1/2 < 0d,2 S 1. This gives the following conditions: 1 2 <0¢I4=s< 1, (a.17) 2 O¢251<=>32§' For example it is found that s = 2/ 3 and 70 = 1 whenever 0 = 04, = 1; therefore 2 (0, yo) = (l, l) is the minimum point. It is found that 2/ 3 < s < l and 70 < 1 whenever 1/2 < 04,2 < l ; therefore (a. 16) gives a minimum point which is not the end point (0, yo) = (1,1) (see Figure 10). Thus, if 0 5s 3 2/3, then there doesn’t exist an internal mini- mum point in the (0, yo) -plane on the curves of constant s; if, however, 2/ 3 < s < 1, then there exists an internal minimum point. A.2 The hard device energy function A12 for material A From (25) and (a. 10)2 the hard device energy function A12 is given by l 1 02 1‘,-,(0) 70+(20-l+1n(20-1)) 'Yo-IXO) A12 = A1203, 7) — 2 I‘,,(0)- l‘,(0)+ 2 1‘”(0) - F,(0)° (“8) Taking derivatives with respect to 70 along the curves of constant 0 gives 8A12 — > 0. (a.19) 870 0 This indicates that the hard device energy function A12 monotonically increases along the curves of constant 0 as 70 increasing from 0 to 1/ (20 - 1) . It is also found that 161 8A — > 0 when 70 2 l. (a.20) This implies that hard device energy function A12 monotonically increases on the line of constant 70, yo 2 l, as stress increases from 1/2 to 1/ 2 + 1/ (270) . Moreover the hard device energy function A12 increases from 1/ 2 to 60 (yo) and decreases from 6% (ya) to 2 2 o on the line of constant 70, with 70 < 1. Note from (48) and (a.4) that the curve $2 is given by 40-1 1 ,_ 0 whenever 1/ 2 < o < 1. From (123) and (a.24) the first derivative of ink!) is given by l 20—1<0’ 1fl/21. Note also that the constant CC = 1/2 and the behavior of strain-softening 611(7) —9 1/ 2 when 7 -) oo which are similar to that of material A. The first derivative of 6” (y) is given by 6,,'(y) = L1>0, (b.2) 73 which shows the different strainosoftening from that of material A. The second derivative of 611(7) is given by - n 3 O” (Y) = _>O» (b3) 7‘ which shows the different curvature restriction. From (b. 1)2 the inverse function of 611(7) is given by I‘ ,(o) = ———1-—— (b4) ’ (20 - 1) V2 From (7) and (b.l)2 the energy density function W (7) for 7> 1 is given by _ 1 ‘Y _ i . W(y) — 2+ 2 27’ 1fy> l. (b.5) Thus from (b4) and (b.5) the energy density function W (1""(o)) can be rewrittenas 163 164 W”(I‘”(0)) = 3(1+j2;__1-./2o—1), if I‘”(0)>1. (b.6) One of the main differences between this material and material A involves the limit defined in (10). Note from (10), (b4) and (b.5) that the limit of WC for material B is given by _ . _ . 1_1 _1 Wc-YlgnJWW) vac) - .19.: '27 - 5 (bi) while recall from (a.7) that the limit for material A is given by We = 00. Hence the area under the strain-stress curve that is also above the horizontal line 0 = 0 c is finite. Accord- ing to (b4), the equation (23) which involves the boundary conditions u(O) = 0 and u(L) = 80 can be rewritten as 0s+0(l --s) = 0 = 70, (1,1)-case, 1 03+ (1 —s) = , (1,2)—case, J20 - 1 Y” 1 (b.8) s + 0 (1 — s) = , (2,1)-case, J20 — 1 Yo I l 1 s + (1 - s) = = , (2,2)-case. J20- 1 J20- 1 J20— 1 70 From (b.8) it is found that . Yo ‘ 0 312(0’ 70) = 1 9 - 0 20—1 7 _ 1 (b.9) . 0 J20 — l s21(0.70) = 1 . — J20 - 1 while 311(0, 70) and 322(0, 70) do not have meaning too. Alternatively (b.9) can be found directly from (25) and (b.4) using F,(0) = 0. From (24), (b.4) and (b6) the hard device energy function A for material B is given by 165 A11 = :2 = A (0, 70) (1,1)-case, 2 11 02 1+ ,__2(:_l-J20-1 .. A12 = —2-s+ 2 (1 —s) = A12(0,70) (1,2)-case, 1+ 2C:_1-./2o--1 . (b'lo) A21 = 2 5+ —2- (1 -s) = A21(0,yo) (2,1)-case, 1+ 1 — 20-1 A22 = 20- i = A22(0,'yo) (2,2)-case, where s is given by (b.9). The soft device energy function E for material B can be obtained from (26), and (b. 10) E.“ = A“(0,'yo)-0yo, (1,1)-case, a = A (0,7 ’07.» (1,2)-case, 12 A 12 a) (b.11) 5.21 = A21(0,70)-0yo, (2,1)-case, :52, = 822(0, 701—070. (2,2)-case. B.l The solution region {212 for material B The material properties and the energy functions are different between material A and ma— terial B. How about the solution regions? From (29) and (b.l) the solution region for mate- rial B is given by 1 . 1 1 1 . £212: {(0,yo)|§<0 smm the location of the minimum point on the curve of constant s can be found by (l-S) -————————— = 0. (b.13) S (20-1)” 4 Thus the root of (b.13) is given by 166 - _ 1 (1 _ s) 2/3) 0 — 0¢2(s) _ §(1+ ( s ) . (b.14) Hence from (47), (b.4) and (b. 14) the minimum 70 on the constant s curve is found to be s 3(s)1/3(1_S)2/3 Yols = 5 2 57¢st (b.15) Therefore, from (b. 14) and (b. 15) for a curve of constant s > 3mm, the minimum point is given by (0.7,) = (04,2(3). 1,261) = (1+1 “-s))” s+3(S)1/3(l—s)2/3) (b.16) 2 2 ’§ 2 . Note that the stress in 912 must obey % < 003 s 1. This gives correspondingly l<0 ¢=ss__l_. oz 2 As an example it is found that s = 1/2 and 70 = 1 whenever 0 = 0(1)2 = 1; therefore (0, 70) = (l, 1) is the minimum point. It is found that 1/2 < s < 1 and 70 <1 whenever I/ 2 < 002 < 1 ; therefore equation (b. 16) gives a minimum point which is not the end point (0, yo) = ( 1, 1) (see Figure 10). Thus, if 0 .<. s 3,1/2, then there doesn’t exist an internal minimum point in the (0, 70) -plane on the curves of constant s; if, however, 1/ 2 < s < 1, then there exists an internal minimum point. The above behavior is similar to the behavior of material A. 167 B.2 The hard device energy function A12 for material B From (25) and (b. 10)2, the hard device energy function A12 is obtained by 1 (1+———J20-1 -I‘ 0 A _ A (0 0) _ 02(I‘”(0)-'Yo) + J20-l )(YO ’( )) (b18) 12 1?- ” 20,,(01- r,(o)) 2 (rum) - 13(0)) ' Taking derivatives with respect to 70 along the curves of constant 0 obtains 8A ’2 > 0. (b.19) 0 0 This means that the hard device energy function A12 monotonically increases along the curves of constant 0 as 70 increases from 0 to 1/J20 - 1. It is also found that 8A12 E Y > 0 when “ya 2 1. 03-20) This implies that the hard device energy function A12 monotonically increases form 0 c to 1/ 2 + 1/ (27:) on the line of constant 70 2 1. However the hard device energy function A12 increases from 1/ 2 to 69 (70) and decreases from 6Q (70) to 0 on the line of constant 70 < 1. 2 2 Note that from (48) and (b.4) it is found that <0 0 whenever I/ 2 < 0 < 1. From (123) and (b.24) the first derivative of f12(0) is given by 1 f12'(0) = O- ‘55—: which means that the driving traction [32(0) decreases as a function of 0 whenever <0, l/2<0<1, (b.25) 1/2 < 0 < 1. Note that the rate of change of the driving traction [32(0) for material B is lower than that for material A. Appendix-C: Ordering of the solution paths in Figure 33 In this Appendix we show that a“ > ab = 7(0, or“) > 7(0, (1") at each fixed 0 in oc < 0 < 1 under the loading given in (149): yo = 1+ at. (c.1) From (142) the loading 70(t) is related to 0(t) and s(t) by 70(1) = 0(t) + T(U(t)) (1 - 5(0). (62) Taking the partial derivative of 1(0, (1) with respect to the loading rate a from (c.2) with fixed 0 gives 0 35 (C3) 3 = «0(0):; 0 0 Similarly taking partial derivative with respect to the loading rate a from (c.1) with fixed 0 gives 8—10- = H a-ai . (c.4) a 0 Ba 0 Recall from (148) that s is given by s = 3(0) + kjgf12(c(t))dt, thus 8s - at 33 o = kf12(0(t))5a a. (a.5) Now (c.4) can be rewritten, with the help of (C5), as a as En, 35 ., T = + -—.——— . (C.6) 0‘ 0 kf12(0(t)) From (C3) and (C6) it is found that 37 ('41) Ta 0 33 = + - . (c.7) 0‘ 0 kf12(o(t))T(o(t)) 169 170 . 9Y0 Solvrng (c.7) for — gives Ba 0 3Y0 I SE a = a . (0.8) + - kf12(o(t))T(o(t)) If we consider t> 0 then --—°— > O. (c.9) act a We note also that (C3) and (C5) now also give that a at 31 <0, _ <0, (for t>0). (c.10) a 0 Ba 0 From (c.9) we can conclude that yo >70 if a">a”. (c.11) a“ or" O O Hence the solution path with loading rate a“ is always above the solution path with loading rate ab if a“ > on”, except at the starting point (0, Yo) =(1,1) where t = 0. -' " “W Appendix-D: Local analysis of the solution paths in Figure 33 near 1:0 for material A. In this Appendix we show the local analysis of the solution paths of Appendix-C near t=0 for Material A. From (139) we can obtain (T(0(t)) + (700) - 0(0) T'(6(t))) (30) -k (0p,(0(t)) - 0(1)) 7310(0) = Y0(t)T(O(t))-(d.1) Consider the loading given in (149): 70 = 1 + at. (d.2) From (a.4) we have that T(0(t)), for Material A, is given by T(0(t)) = Fume» - om = 3373—3 - om. (d3) Since 0(t) = 1 when t is near zero, we let p=sm-1 we so that p = 0(t). (d.5) Then T(0(t)) in (d.3) can be expanded to give mm» = To» = - 3p + 0022). (do and flow) = 73(1)) = 9p2 + 0(p3). (d7) The first derivative of T(0(r)) is given by T’(c(t)) = T'(p) = — 3 + 0(p). (d.8) From (119), (121) and (130) we have fume» = (o,,(o(r)) - am) How) = ;(,,T(v)dv. (a.9) Therefore with the help of (d.4) and (d.6) we have (o,,(o(t)) - om) mm» = 137(va = - 319 + 0 (p3). - «1.10) Now (d.1), with the help of (d.2), (d.4), (d.5), (d.6), (d.7), (d.8) and (d.10), can be rewritten 171 172 as 1—3p+0(p’)+ (at-p) (—3+0(p))1p 3 (d.11) = Ice-3125009)) [9p2+0(p3)1+a(—3p+0(p2». Since lpl « 1, if a it 0 then (d.11) is dominated by -3atp’ = -3ap. (d.12) Differential equation (d.12) for p(t) is to be solved subject to p(0)= 0(0)-1=l-1=0. There- fore the leading order solution for p(t) is ' P0) = C ° t. (d.13) where c is an undetermined constant. Hence from (d.5) we obtain 0(t) = c + 0(1), 0(t) = 1 + ct+ o(t). (d.14) Therefore the initial slope of the solution path at t = O is dy_or E. 35 — (d.15) The value of which is not provided by the present analysis. This slope is, however, related to the constant s curves in 912 which the solution path follows as it leaves (00, yo) =(1,l), since (25) now gives 1 (l s(t)=§(2+F)+0(l). (d.16) Thus c is given in terms of the loading rate a and the initial nucleation site so, by (I (d.17) andwenotewith a>OthatC>Oif2/3oo and 50—) - Strans’ WIN If we consider an extended expansion p(t) = ct+dt2 +ez3 + gt4+ 0 (t4), (d.18) 173 then (d.1) can also be expanded to yield (- 48 + 4otc2 — 3da) :2 + (16c4— 16c3ot — 16c2d+ 16cda - oea) :3 + (— 48:? + 48oto4 + 13.51tc4 + 80dc3 — 72dotc2 - 20cc2 + 24eotc - 20ch (d.19) - 9af+ 12(12(1)t4 + o(t“) = 0. According to the coefficient of t2 in (d.19) we have 43 a _4a2 3 0‘ (3s0-2)3 (d.20) Thus the coefficient d is dependent on the value a and so, but is independent of k. Accord- ing to the coefficient of t3 in (d. 19) it is found that c3 (4c2-soc+o2) _ 80960-1) (so-2) 2 8 e:— 9 a (3s0-2)5 (d.21) Hence the coefficient e is also dependent on the value a and so, but is independent of k. According to the coefficient of t4 in (d.19) we have -1 c‘ (640c3 - 96Oac2 + 288a2c - iiiltot2 + 32013) g = 371' 3 = (l 3 (d.22) 7(32(1(s0-1)(.33‘24-230-4)-+-3k(-27S:+5453-3630-1- 8)). 2(3s0-2) Therefore the coefficient g is dependent on the value a, so and k. Hence from (d.5) we obtain 0(t) = c + 2dr + 3c:2 + 4gt3 + o(t3), 0(t) = 1 + ct+ d:2 + et3 + gt4 + o(t“), (d.23) where d, e, and g are given in (d.20), (d.21), and (d.22), respectively. Hence from (25), (d.2), (d.23), (d.6) we can find _1 or 1 2 _ 1 [“12 3 3 s(t) _ 5(2+?) -§kc t3+o(P) _ SO-EWt +0(t ). (d.24) Note with the help of (d.9) and (d. 10) that (d.24) also can be found by integrating s = kf. List of References List of References [195613] J.D. Eshelby, The continuum theory of lattice defects, Solid State Physics (Edited by F. Seitz and D. Turnbull) vol 3 , Academic Press, New York, 79-144. [197OE] J.D. Eshelby, Energy relations and the energy-momentum tensor in continuum mechanics, Inelastic Behavior of Solids (Edited by M.F. Kanninen et al.), MacGraw-I-Iill, New York, 77-115. [1973D] C.M. Dafermos, The entropy rate admissibility criterion for solutions of hyperbol- ic conservation laws, J. Differ. Equations 14, 202-212. [1973L] P.D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves, CBMS Regional Conference Series in Applied Mathematics, 1 I , SIAM Publications: Philadelphia. [1975B] J.L. Ericksen, Equilibrium of bars, J. Elasticity 5, 191-201. [1975513] J.D. Eshelby, The elastic energy-momentum tensor, J. Elasticity 5, 321-335. [1975R] J .R. Rice, Continuum mechanics and thermodynamics of plasticity in relations to microscale deformation mechanisms, Constitutive Equations in Plasticity (Edit- ed by AS. Argon), MIT press, Cambridge, MA. 23-79. [1977K] J.K. Knowles, The finite anti-plane shear field near the tip of a crack for a class of incompressible elastic solids, International J. of Fracture 13, 611-639. [1978K] J.K. Knowles and E. Stemberg, On the failure of ellipticity and the emergence of discontinuous deformation gradients in plane finite elastostatics, J. Elasticity 8, 329-379. [19798] Z.P. Bazant and 8.8. Kim, Plastic-fracturing theory for concrete, J. Eng. Mech. Div., ASCE III, 407-428. [1979K] J.K. Knowles, On the dissipation associated with equilibrium shocks in finite elas- ticity, J. Elasticity 9, 131-158. [1980J] R.D. James, The propagation of phase boundaries in elastic bars, Arch. Ration. Mech. & Anal. 73. 125-158. ‘ [19800] D.R.J. Owen and E. Hinton, Finite elements in plasticity: theory and practice, Pin- 174 175 eridge Press Limited, 13-31. [1982K] N. Kikuchi and N. Triantafyllidis, On a certain class of elastic materials with non- elliptic energy densities, Quarterly of Applied Mathematics 40, 241-248. [1983A] R. Abeyaratne. An admissibility condition for equilibrium shocks in finite elastic- ity, J. Elasticity 13, 175-184. [1983H] R. Hagan and M. Slemrod, The viscosity-capillarity criterion for shocks and phase transitions, Arch. Ration. Mech. & Anal. 83, 333-361. [1984R] H.E. Read and G.A. Hegemier, Strain softening of rock, soil and concrete - a re- view article, Mech. of Material 3, 271-294. [19858] Z.P. Bazant and TB. Belytschko, Wave propagation in a strain-softening bar: ex- act solution, J. Eng. Mech. Div. ASCE III, 381-389. [19850] M. Ortiz, A constitutive theory for the inelastic behavior of concrete, Mech. of Ma- terial 4, 67-93. [1986H] H. Hattori, The Riemann problem for a van der Waals fluid with entropy rate ad- missibility criterion - isothermal case, Arch. Ration. Mech. & Anal. 92, 247-263. [1986?] TJ. Pence, On the emergence and propagation of a phase boundary in an elastic bar with a suddenly applied end load, 1. Elasticity I 6, 3-42. [1987A] R. Abeyaratne and J .K. Knowles, Non-elliptic elastic materials and the modeling of dissipative mechanical behavior: an example, J. Elasticity 18, 227-278. [1987AA] R. Abeyaratne and J.K. Knowles, Non-elliptic elastic materials and the model- ing of elastic-plastic behavior for finite deformation, J. Mechanics and Physics of Solids 35 . 343-365. [1987B] T. Belytschko, XJ. Wang, Z.P. Bazant and Y. Hyun, Transient solitons for one- dimensional problems with strain softening, J. Applied Mechanics 54(3), 513- 518. [1987?] R.L. Pego, Phase transitions in one-dimensional nonlinear viscoelasticity: admis- sibility and stability, Arch. Ration. Mech. & Anal. 97, 353-394. [1988A] R. Abeyaratne and J. K. Knowles, On the dissipative response due to discontinu- ous strains in bars of unstable elastic material, International J. of Solids and Structures 24, 1021-1044. ' [198 8AA] R. Abeyaratne and J. K. Knowles, Unstable elastic materials and the viscoelastic 176 response of bars in tension, J. Applied Mechanics 55, 491-492. [19888] 8.A. Silling, Numerical studies of loss of ellipticity near singularities in an elastic materials, J. Elasticity 19, 213-239. [198888] S.A. Silling, Finite-difference modeling of phase changes and localization in elasticity, Comput. Methods Appl. Mech. Eng. 70, 251-273. [1989A] R. Abeyaratne and CH. Jiang, Dilatationally nonlinear elastic materials-I, some theory, International J. of Solids and Structures 25, 1201-1219. [1990A] R. Abeyaratne and J.K. Knowles, On the driving traction acting on a surface of strain discontinuity in a continuum, J. Mechanics and Physics of Solids 38, 345- 360. [1991A] R Abeyaratne and J .K. Knowles, Kinetic relations and the propagation of phase boundaries in solids, Arch. Ration. Mech. & Anal. 114, 119-154. [1991AA] R. Abeyaratne and J .K. Knowles, Implications of viscosity and strain gradient effects for the kinetics of propagating phase boundaries in solids, SIAM J. Ap- plied Math. 51, 1205-1221. [1991AAA] M. Affouf and RE. Caflisch, A Numerical study of Riemann problem solu- tions and stability for a system of viscous conservation laws of mixed type, SIAM J. Applied Math. 51, 605-634. [1991B] J .M. Ball, P.J. Holmes, R.D. James, R.L. ?ego, and P.J. Swart, On the dynamics of fine structure, J. Nonlinear Sci. 1 , 17-70. [19916] MB. Gurtin, On therrnomechanical laws for the motion of a phase interface, J. Ap- plied Math. and Physics 42, 369-388. [1991?] T.J. Pence, On the encounter of an acoustic shear pulse with a phase boundary in an elastic material: reflection and transmission behavior, J. Elasticity 25 , 31-74. [1991??] T.J. Pence, On the encounter of an acoustic shear pulse with a phase boundary in an elastic material: energy and dissipation, J. Elasticity 26, 95-146. [1992A] R. Abeyaratne and J.K. Knowles, On the propagation of maximally dissipative phase boundaries in solids, Quarterly of Applied Mathematics 50, 149-172. [1992L] J. Lin and T.J. Pence, Energy dissipation in an elastic materials containing a mo- bile phase boundary subjected to concurrent dynamic pulses, Transaction of the Ninth Army Conference on Applied Math. and Comput., (Edited by F. Dressel), 437-450. 177 [1993L] J. Lin and T.J. Pence, Kinetically driven elastic phase boundary motion activated by concurrent dynamic pulses, Transaction of the Tenth Army Conference on Applied Math. and Comput., (Edited by F. Dressel), 69-86. [1993LL] J. Lin and TJ. Pence, On the dissipation due to wave ringing in non-elliptic elas- tic materials, To appear in J. Nonlinear Science. [19938] S.?. Shah, C. Ouyang, and DA. Lange, Fracture behavior of cement-based mate- rials. Materials Research Society Bulletin, 55-59.