Home

" 1'r_é§O"LR-758 82000 User Manual for Crack

image

Contents

1. 2 8 While the global field u extends itself over the whole domain the singular field u is nonzero only in a small domain around the crack tip The singular crack tip shape functions u are constructed of the classical crack tip solutions derived in the linear theory of elasticity This solution is not directly applicable to the nonlinear case because we have to take into account the rotation of the crack tip in space during the deformation process If the shape functions corresponding to the classical crack tip solutions are given by v the function u gt is given by u R r v 8 0 k R Sp vo 2 9 where R Sr denotes the rotation at the crack tip the location of the crack tip Note that the rotation matrix is dependent on the solution u so that one could write u R p G v 6 0 k OR Sr S vo 2 10 The configuration in the deformed state can be seen to be attained in two stages the first stage by the displacement u the second by the addition of u In this way we identify an intermediate state B described by x X u y amp 0 U amp 2 11 The base vectors that belong to the undeformed state are given by d A y E 0 2 12 i xe 4 14 94 B2000 User Manual page 11 In the intermediate state the base vectors are b x 6 A Sy amp 6 Aj u 2 13 t oe I Note that a comma followed by an index i denotes partial diff
2. where NTE are the finite element shape functions Returning to equation 2 24 we decompose this set of nonlinear equations in the form q Q q 1 2 Figure 4 Separation of the analytical displacement solution corresponding to the stress intensity factors and the underlying finite element displacement field 4 14 94 B2000 User Manual page 15 f g f q k A 0 2 27a h k A O 2 27b where f f e Ry h e R The first term in equation 2 27a corresponds to the contributions of the finite element discretization scheme The two other terms in these equations arise when we apply the enrichment scheme described in the previous sections Equations 2 27a b can be seen to be derived from the requirement that the total potential energy of the structure with its loads must be stationary The contributions f f and h are then computed from oP O OP where the potential energy is given in two parts Protat P q 2 P q k A 2 29 The first term in equation 2 29 refers to the global part of the formulation while the second term P is due to the mode enrichment process 2 1 1 5 Solution of the corrector equations The nonlinear equations 2 27a b are solved for a range of values of the loading depending on the particular requirements for the analysis of the problem at hand This solution process is carried out step by step by variation of the load parameter or more generally by variation of a generalized
3. thus 3 5 3 6 3 7 3 8 3 9a 3 9b 3 9c 3 10 4 14 94 B2000 User Manual page 39 dE s 0 3 11 or Bo f 3 12 where B and f were already defined in equations 3 9a b The continuous stress field o is now obtained by solving equation 3 12 This way of working can be applied for all stress components Next the energy norm per element can be calculated as 1 2 leel Jeg De dv 3 13 y e For the total energy norm this yields 2 lel Fjes 3 14 e Although the absolute value of the energy norm has little physical meaning the relative percentage error is more easily interpreted n Hel x 100 3 15 where 1 2 IUI fo D odv 3 16 v and represents the square root of two times the total strain energy Equation 3 15 is the so called weighted RMS root mean square percentage error in the stresses The percentage error T is calculated for the whole domain and for all elements separately The energy norm for the separate elements gives the contribution of each element to the total error and local mesh refinement is performed on the basis of these element errors A very common requirement is to specify a maximum percentage error for the total energy norm One thus requires that after the final analysis is completed the condition s w a I UL UL 1 A E A T Bh E ln page 40 B2000 User Manual 4 14 94 T lt 3 17 is satisfied for the whole d
4. 1 0 1 0 equal to the closest boundary of this interval If for instance the final values of and gt are 0 5 and 1 05 respectively than gt is set equal to 1 0 Figure 17b 4 1 2 The mapping process The mapping process for the deflections and the base vectors of the orthogonal triads describing the rotations in the element nodes from the old finite element mesh onto the new one is straightforward once the position of the new nodes is determined with respect to the old mesh The coordinates of a projection point P were obtained as x y z N E gt x y z 4 1 page 48 B2000 User Manual 4 14 94 where x y Z p are the coordinates of the projection point P are the parametric coordinates of P P P N denotes the corresponding shape functions and x y z are the nodal coordinates of the element of the old mesh in which P is situated The deflections and orthogonal triads in this projection point P and thus in node n of the new mesh can be determined now in the same way yielding u n N iy u n 4 2 where u n p are the deflections and triads corresponding to the projection point P and u n are the deflections and triads of the nodes of the old element that contains P As already mentioned in the introduction the solution data that are mapped do not correspond to an equilibrium state for the new finite element mesh However when the difference between the new and the old
5. 94 14 15 16 17 18 19 20 82000 User Manual page 31 stress intensity factors and the corresponding integration order INTM 4 only used if MENRI 1 TIELENM 1 NCRITM identification numbers NCRITM I5 of the enriched elements used for the calculation of the membrane stress intensity factors 4 only used if MENRI 1 NCRITBT INTBT number NCRITBT of enriched 2 I5 elements used for the calculation of the bending and tearing stress intensity factors and the corresponding integration order INTBT 4 only used if MENRI 1 and NB1 1 or NB2 1 or NK3 1 IELENBT 1 NCRITBT identification numbers NCRITBT I5 of the enriched elements used for the calculation of the bending and tearing stress intensity factors 4 only used if MENRI 1 and NB1 1 or NB2 1 or NK3 1 Attention Elements that are enriched for the calculation of the bending or tearing stress intensity factors or the elastic T term also have to be enriched for calculating the membrane stress intensity factors see input card 10 NPBRDM number of corner nodes bordering I5 the enriched element domain used for the calculation of the membrane stress intensity factors 4 only used if MENRI 1 IPBRDM 1 NPBRDM identification numbers NPBRDM I5 of the nodes bordering the enriched element domain that is used for the calculation of the membrane stress intensity factors 4 only used if MENRI 1 NPBRDBT number of corner n
6. JFIM sNOTLOIS Z IJIM O a Z ayn OSLIUM PUT INOIQNS pua soro z soro 1 soTo nurT3uoo OI 0306 zypuse sseux sainqjoniys jo SSEW T307 C 3T20 sseux samnjoniys jo ssew e209 4 9 JFIM 303 squawata e113 Jo pue penb jo eaze 18307 a x E O37IM 7032 S3U W T 8113 zo pue penb jo eaze 183707 9 IFIM ano 03 U JJJIM SquawaTe qwnd JUNOD y E IIJIM ano OF UaIITIM SsJU W T OVNO IUNOD y 9 IFIM Z JLINM TTE u y3 66 b l4qul1 J ST ssewx SSEU p 31T1IM wa awx Sseux Ss eux 54 uy 346J m 18703 waTewx weTse LHOTGM TIEI 2038 73038 4 3971s u r e 4o3e 3o3e water 64qur z KA x wadw TIE 00lI 6quz 2 A X 13 11808 1189 Jipu g g 53ur p bU tE Bauz Z 53UF T 53UF 1 pear T pear T qunod qunoo u ua b aur JI JYrpu tp Baug c baup 2 53uy 1 bauy T peas T peas T qunod qyunoo u u3 p bare quy JI uaya z b tu JT asta e t pear eg6 9TaE 38w10 Zquy 2 Z307 A ZT x ST T peas u u3 U b r3ur JI Ey Zquyz Tuy oOT pua T peas OdLlum Tres Q SSPWX On JoO1e O quNeS JIANG IVILNGNOGS SS322V AILLYWYOJ WAO S M3N SnLVLIS I3WVN 3TI3 TLNOUN LINN NIJO 3513 TLNOUN GNIMIY IVILNGNDGS SS322V GALLVANOd W3O4 qTO S 01VLS TANWN 31T 14 TIQOSN LINQ0 N3dO NHL ISIX33 dI IS
7. Vol 28 1989 pp 1461 1473 A P KFOURI Some evaluations of the elastic T term using Eshelby s method Int J of Fract 1986 pp 301 315 f H A J KNOPS Computation of stress intensity factors by the mode enrichment technique with applications to geometrically nonlinear problems Report LR 736 TU Delft Faculty of Aerospace Engineering Delft The Netherlands 1993 H A J KNOPS Numerical simulation of crack growth in pressurized fuselages Ph D Thesis TU Delft Faculty of Aerospace Engineering Delft The Netherlands 1994 H A J KNOPS Input description for the QUADMESH mesh generator Report LR 759 TU Delft Faculty of Aerospace Engineering Delft The Netherlands 1994 S MERAZZI B2000 a modular finite element analysis system A summary NLR MEMO SC 85 041 U National Aerospace Laboratory NLR The Netherlands November 1985 S MERAZZI B2000 Data access methods and database description handbook version 1 2 SMR Corporation P O Box 41 CH 2500 Bienne 4 Switzerland September 1985 C D MOTE Jr Global local finite elements Int J for Numer Meth in Engrg Vol 3 1971 pp 565 574 D R J OWEN and A J FAWKES Engineering fracture mechanics Numerical methods and applications Pineridge Press Ltd Swansea United Kingdom ISBN 0 906674 26 3 1983 E RIKS Some computational aspects of the stability analysis of nonlinear structures Comp Meth in Appl Mech and Engmg Vol 47
8. crack tip will give a considerable erroneous contribution to the total in plane strains This error can be corrected for by subtracting the analytical in plane strains from the total strains before calculating the stresses see section 2 1 1 3 However in order to get hold of these analytical in plane strains the orientation of the crack tip zone after rotation has to be determined For this purpose the elements IEROT 1 2 and nodes IPROT 1 3 are used These elements and nodes have to be defined as pictured in Figure 12 IEROT 1 TEROT 1 IPROT 2 IPROT 2 IPROT 1 IPROT 1 IPROT 3 IEROT 2 IEROT 2 4 node Bathe element 8 or 9 node Bathe element Figure 12 Definition of the elements IEROT 1 2 and the nodes IPROT 1 3 6 The nodes and elements along the crack faces are used to calculate the elastic T term via a least square fit technique see section 2 1 2 The nodes NODEUP 1 NUMNOD and NODELW 1 NUMNOD and the elements IELEUP 1 NUMNOD and TELELW 1 NUMNOD are defined as pictured in Figure 13 2 3 2 Input file config Information concerning the configuration number that the separate elements belong to see point 1 and 2 in Section 2 3 1 is given in input file config The 2 input cards that are needed in this input file are as follows card description format 1 NFIG number of element configurations involved 2 I5 in the total mesh NELE total numer of elements 2 IELE elem
9. mode enrichment In section 2 1 1 3 it was shown that the difference between the intermediate state B and the final deformed state B is furnished by the inclusion of the analytical displacement solution u in the total displacement field This analytical displacement solution is only locally defined and therefore its influence on the total displacement solution is mainly restricted to a small area around the crack tip Actually the enrichment process represents a linear technique applied to a small part of a geometrically nonlinear deformed structure which gives rise only to small changes in the general solution Therefore one can apply mode enrichment as an additional step after the nonlinear state of deformation of the structure described by the global field u is completed It turns out that because of the linear character of the mode enrichment process the stress intensity factors and thus the final deformed state is obtained within a few iterations only This technique is particularly useful in cases where the knowledge of the stress 4 14 94 B2000 User Manual page 19 intensity factors is not required at each point of the load deformation path but only for the final load value There is yet another reason why this a posteriori technique is attractive to use As will be described in section 2 1 1 7 a higher order Gaussian integration scheme has to be applied for the enriched elements to get accurate solutions for the stress intensity factors Thi
10. order of this numerical integration procedure has to be chosen in such a way to achieve an acceptable level of accuracy for the finite element results including the stress intensity factors A widely used numerical integration technique is the so called Gaussian quadrature procedure 1 In this procedure the position and the weight of the integration points is optimized in such a way that the finite element matrices containing the polynomial interpolation functions for the isoparametric elements are integrated exactly This exact integration can be achieved with a low order 2x2 or 3x3 integration scheme depending on the type of element used 4 or 8 9 nodes The enriched elements however also contain shape functions corresponding to the analytical crack tip solutions for the different modes of cracking Due to these additional shape functions one can not suffice with a low order integration scheme for these elements 2 8 11 12 19 25 26 28 In 14 15 a few example calculations are presented and it is shown that a higher order integration scheme of 10x10 will give converged solutions page 20 B2000 User Manual 4 14 94 2 1 2 The elastic T term via mode enrichment 2 1 2 1 Introduction The elastic T term is calculated by extracting it directly from the finite element displacement field around the crack tip In the literature one can find several other techniques for calculating the elastic T term However these methods make use of addition
11. that one should enforce this condition at the boundary of the enriched area A simple way to accomplish this is to use a transition function This function is equal to 1 0 for the enriched elements next to the crack tip For the transition elements though this function is equal to 1 0 next to the enriched elements and equal to 0 0 next to the conventional elements Figure 2 while the descent of this transition function is assumed to be bilinear How to enter all this information in the B2000 routine b2contcrack will be discussed in section 2 3 1 The particular application of the mode enrichment technique that we have in mind is the problem of a longitudinal crack in a pressurized fuselage This implies that we are enriching a curved surface with analytical crack tip solutions that only hold for a flat crack tip zone However this difficulty is circumvented by restricting mode enrichment to a small area around the crack tip which is approximately flat Therefore the possible curvature of a structure plays no role as far as the mode enrichment process is concerned 2 1 1 3 Corotational formulation The set of singular functions discussed in section 2 1 1 1 is taken from the linear elasticity solution for an infinite plate with a crack under general loading conditions It is defined with respect to the orientation of the undeformed crack geometry which in the linear theory does not appreciably change during the loading phase When the problem is geo
12. the elastic T term 26 NODELW 1 NUMNOD identification numbers NUMNOD I5 of the nodes along the lower crack face that are used for the calculation of the elastic T term 27 TELELW 1 NUMNOD identification numbers NUMNOD I5 of the elements along the lower crack face that are used for the calculation of the elastic T term t 1 The orthogonal triads in the element nodes are calculated in subroutine cornrd First the normal vector is calculated The next step is to calculate 1 of the 2 in plane rotation vectors by determining and normalizing the vector product between this normal vector and 1 of the cartesian unit vectors in cornrd e gt 0 1 0 is used for this purpose f However with this technique one can end up with in plane rotation vectors which are not aligned with the model boundaries and this can pose some difficulties if one wants to apply boundary conditions along these boundaries However the options IFIXV1 1 or IFIXV2 1 offer the opportunity to fix 1 of the 2 in plane vectors on beforehand For a cylindrical coordinate system for instance one can force 1 of the 2 in plane rotation vectors to be aligned with the axis of this cylindrical coordinate system Because the 4 14 94 B2000 User Manual page 33 normal vector in a cylindrical coordinate system is automatically forced to point in the outward radial direction the second in plane rotation vector will be aligned with t
13. v vv E A Eq vk 8 Ak en 2 20 where we assume that the base vectors Ak cp at the crack tip for the undeformed state are unit vectors and mutually orthogonal Substitution in equation 2 19 yields 1 T T G j u Rv 5 tic RAg or Yki Yki RAK cry Pj cer VkiVkj 1 2 21 The base vectors bir are obtained from the base vectors Ay cr by a rotation which we identify with R and a small deformation caused by the field u It is now natural and consistent with our previous assumptions to neglect the deformation of the triad Ag cr and replace bi cr by RAk co The additional strains Gi that arise can then be simplified to the expressions 1 I Gi u u G u Rv 5 Vij F Vii f Vi Vk 2 22 Comparing this additional Green Lagrange strain Gi with the truncated strain field derived on the basis of linearized elasticity theory one can see that the only difference lies in the nonlinear term v vk The effect of this nonlinear term on the total additional strain Gi remains small as long as the displacement gradients are small However when approaching the crack tip these displacement gradients will increase and so will the nonlinear term Vk Vk This increase of v V will cause a discrepancy between the truncated strain field and the Green Lagrange strain field f The concept of singular stress and strain field around the crack tip is derived on the basis of linearized elasticity theory which is formu
14. with the old mesh PBS may not be increased during the new load history PBS the external forces that were originally defined in load case 1 now have to be defined as load case 2 In order to apply the mode enrichment technique the b2contcrack processor needs some special information concerning the enriched element zone As an additional feature one has the opportunity to use a cylindrical coordinate system in which to define the external loads and boundary conditions Of course this implies that the nonlinear analysis is also performed with respect to this cylindrical coordinate system In what follows the separate input cards are described in a brief way For an explanation of the use of some of these cards the reader is referred to the end of this section card description 1 ICRSYS used coordinate system ICRSYS 1 cartesian coordinate system ICRSYS 2 cylindrical coordinate system 2 CYLC 1 3 1 cylinder axis only used if ICRSYS 2 3 CYLC 1 3 2 point on the cylinder axis only used if ICRSYS 2 l 4 IFIXV1 IFIXV2 initial in plane vectors for the definition of the rotations are either defined on the next input card or determined in the element format 15 3 E12 5 3 E12 5 2715 page 30 10 11 12 13 B2000 User Manual 4 14 94 processor only used if ICRSYS 2 IFIXV1 1 rotation vector 1 is defined on the next input card IFIXV1 0 rotation vector 1 is determi
15. x Cot Cix Coxe O 2 46 In this equation Co represents the elastic T term while C4 Cz etc correspond to the higher order terms Further x denotes the distance to the crack tip ell LEAL Aphendscwncve si wt edhe ed Serna bar page 24 B2000 User Manual 4 14 94 Figure 7 Crack tip before and after rotation According to the least square fit technique the error between this continuous function for the stress field and the discrete stresses calculated in the elements behind the crack tip has to be minimized where this error is defined as n E o x 9 2 47 i In this equation x denotes the distance of the midpoint of the element sides along the crack face to the crack tip while o x and o are the values of the continuous stress field 6 x and the discrete finite element stress at point x respectively Figure 8 Further n represents the number of elements over which this least square fit is performed Minimizing the error E error with respect to the coefficients C in the polynomial function yields G xc Kn x x crack tip Figure 8 Extrapolation of the element stresses to the crack tip via the least square fit procedure 4 14 94 B2000 User Manual page 25 0 2 48 for all coefficients C This results in a matrix equation the order of which is determined by the number of coefficients in the continuous stress function o Finally the elastic T term Co in equation 2 46 is
16. 1984 pp 219 259 E RIKS and H A J KNOPS Mode enrichment and the analysis of through cracks in thin walled shells under internal pressure Report LR 736 TU Delft Faculty of Aerospace Engineering Delft The Netherlands 1993 T L SHAM The determination of the elastic T term using higher order weight functions Int J of Fract Vol 48 1991 pp 81 102 G STANLEY I LEVIT B HURLBUT and B STEHLIN Adaptive refinement AR strategies for shell structures Part 1 Preliminary research Lockheed Contract Report LMSC F318483 May 1991 G STRANG and G J EIX An analysis of the finite element method Prentice Hall Inc Englewood Cliffs New Jersey 07632 1973 D M TRACEY Finite elements for the determination of crack tip elastic stress intensity factors Engrng Fract Mech Vol 3 1971 pp 255 265 M L WILLIAMS The bending stress distribution at the base of a stationary crack J of Appl Mech Vol 28 1961 pp 78 82 ka alat ae 4 14 94 B2000 User Manual page 57 28 W K WILSON Combined mode fracture mechanics Ph D Thesis University of Pittsburgh 1969 29 O C ZIENKIEWICZ and J Z ZHU Mesh regeneration and automatic adaptive analysis in Computational Methods in Applied Sciences by Ch Hirsch Elsevier Science Publishers B V ISBN 0 444 89795 X pp 13 24 When ell JUMAH piii iced nal a a bgt TETTA REPRE way owx WoTee LHOIGM urqanoiqns pus ZOU
17. 4 14 94 3 Attention Configuration 1 always represents the shell surface for which the orthogonal x 1 triads calculated as described under point are correct These vectors define the crack tip coordinate system that is used to calculate the stress intensity factors and the elastic T term Figure 10 The strains in these elements are calculated with respect to this crack tip coordinate system before transforming them to the global coordinate system CRADP 1 3 CRADV 1 3 crack tip Figure 10 Definition of the crack tip coordinate system 4 An example for the definition of the enriched elements the bordering nodes and the integration order is given below Figure 11 NCRITM 16 INTM 10 TELENM 123456789 1011 12 13 14 15 16 NPBRDM 17 IPBORDM 123456789 1011 12 13 1415 1617 card 13 NCRITM INTM card 17 NPBRDM card 14 IELENM 1 16 card 18 IPBORDM 1 17 Figure 11 Definition of the input data for the enriched element zone used for the calculation of the membrane stress intensity factors the same type of definition also holds for the bending and tearing stress intensity factors 4 14 94 B2000 User Manual page 35 These corrections are due to the fact that the in plane strains according to the analytical crack tip displacement fields will go to infinity when approaching the crack tip Because of the quadratic terms in the Green Lagrange strains these large analytical strains close to the
18. AZ UTS 6 Zoe ST ZEUUaZ UTSeZaabT 3 EIUJ ZUTS Zv Zt TTUdeZ UTSaZee Il GZ swayoe s r5uerz qa omy Fo jo eaze ayy jo wns ay se adejpe JO ePaze Jo uorqelnoIeo 2 JIpu ST T bT_EdU S098 7 47 TebT pe Tdup sooe z pyd asta O TSYR o ztud 3y3u u lI ieI nbuer13 1037 uor3ad ox 5 uayr 6qup bes p OUP J eT 2T T_ dut sose ay T TI EZ Tduy sooe pyd 2 p2 2 1Z 3 g4 pA s CA TA _ 9 x bx x TX bI gedu 2 22 x Z 12 3 g amp z 0 gA 1A _ 3 x x Ex Tlx Zl eduy TZ pZ 1Z 2 3 TA pA x TA CA _ 9 Ix px TX x pe dur 12 2 12 22 3 T E T Z _ 3 TX X x Ix zx z duy s jbue pasojToua jo uoT3eInoIeo fo Go us Zex 2 b2 Z s K PKA Z x pX SGT S as s 12 2 Z IK bpA Z s IX bpX PT Gow Z Z 2 Za TA A Z TX EX HET Gove Sun Z2 2 Z uu ZA EA 2 ww ZK EX Z1 S as Zee 12 22 24 0 TA Z4 2 ue TK ZX TT quoweTapenb jo euoberp pue s pts jo suq5u rT Jo uot zepnorteos 5 pOQUT Z p2 poquT poqUyT X px g 343u1 z Z go3ur A P pu g523ur1 x g x u r e u Tr e umn u l ux z53ul1 zZ zz 5y qJusuwete jo sseu 2 z523u1 A Z amp JWETe wx 31092 WxeWaTee ux ZOqUyT x zx 0602 FWETOUT 06 0U4S Sb Ix TWE TOUTS pS pusyweTe ux T23 uy z Tz 31027 1I0JOQ41 102V ux T22UF 1 i Z w
19. E Repon LR 758 B2000 User Manual for Crack ta Analysis Solution Data Mapping and Error Calculation April 1994 H A J Knops ihr T U De ift Faculty of Aerospace Engineering Delft University of Technology B2000 User Manual for Crack Analysis Solution Data Mapping and Error Calculation H A J Knops 4 14 94 82000 User Manual Contents 1 Introduction 2 Crack analysis 2 1 Theoretical background 2 1 1 Mode enrichment 2 1 1 1 Basic principle 2 1 1 2 The enriched element domain 2 1 1 3 Corotational formulation 2 1 1 4 The resulting nonlinear equations 2 1 1 5 Solution of the corrector equations 2 1 1 6 A posteriori mode enrichment 2 1 1 7 Numerical integration 2 1 2 The elastic T term via mode enrichment 2 1 2 1 Introduction 2 1 2 2 Solution strategy 2 1 2 3 Numerical procedure 2 2 Interactive input commands 2 3 Additional in and output files 2 3 1 Input file mode_enrich 2 3 2 Input file config 3 Error calculation 3 1 Theoretical background 3 1 1 The energy norm 3 1 2 Modified version of the energy norm 3 2 Interactive input commands 3 3 Additional in and output files page 1 page oN Q Q Q wD 15 18 19 20 20 21 22 25 29 29 35 37 37 37 41 42 page 2 B2000 User Manual 4 6 3 3 1 Output file ERRORS Solution data mapping 4 1 Theoretical background 4 1 1 Node projection 4 1 2 The mapping process 4 2 Interactive input commands 4 33 Additional in and output files 4 3 1 Input
20. IXJJ LSIX TIWYN JTI d3 UINONT ano d o3ji3 aueu g l3nozu UMOU XUN snae3s qno lI1jJ 31un u do pro snae3s ur I1 T 3 Jun u do ano 0zZ89 Pear a r1jandono Jo aweu z 3u3g s JJIM UT a OZB a e pear a saTFF ynduy Jo aweu 134Ug e 93720 spxazy IeOTDOI we Taux sseux water yo e 0000T 2 00001 A 00001 x TeAaZ T3anozuqunoS bp 63u1 3uy aur Ty rabaquy OZ Il ueu0z 3no OZ U 1937981842 PUVUUSCR UCT UCU TTSCTIC VOCS ESOC ELE LESS SES SS eee 66T OT 62 uasuei3 S ZLGWOW UT pasn aq ued YOTYM ATTJ jxq e UFeIGO o3 pasn aq uayy ued q z Terjneu p ljrpou Huy j Nsez ql payetep aq oq aaey 66 IO Z T uey STaQUTOd 1 470 YITH 6buyuuybaq spzed TIE STTI STUI UI NWHLVd WOTJ STF Z Hd TerINeu e SI aTTI qndut PUL ZTQWDV ot NWULvd uo3r3j werhord 5ej3 qur ue S GWOLvd WET UPCTOPVECTEUSCCOCLOOSOSCECCCOCEL ESET SESE eee z o y e Tear BPO TTauy aqwDiWa uei5o zd Sl ot ovuvuVv oo UoOD0 UY lili
21. ack tip solution in a geometrically nonlinear environment with large rotations and translations Finally because the governing equations are nonlinear and because it is advantageous to keep the required changes in the solution algorithm at a minimum it is necessary to discuss the way the governing equations are solved in the altered situation 2 1 1 2 The enriched element domain The first aspect which needs attention is the size of the domain in which mode enrichment is applied As discussed in 14 15 the crack tip solutions for the two bending modes and the tearing mode are only valid within a small boundary layer Therefore in contrast to the page 8 B2000 User Manual 4 14 94 Z ZZZ ZZZ L ENNY ISL _ LZ LLY Z lt A 1 0 Z S _ Z NS sosssssssssssssssss 0 0 crack tip Figure 2 Definition of the bilinear transition function membrane stress intensity factors mode enrichment for the bending en tearing modes has to be restricted to a small boundary layer in order to obtain accurate results for these stress intensity factors This is the main reason why the enriched element domain for the membrane modes on one hand and the bending and tearing modes on the other hand can be chosen independently Furthermore we will restrict the mode enrichment process to some small area around the crack tip In a strict sense one should require that the summation of u and u remains Co continuous everywhere in the domain In practice this means
22. al load systems 7 13 or weight functions 23 and therefore involve much more effort to determine the elastic T term especially when dealing with a geometrically nonlinear problem The analytical plane stress displacement solution around the crack tip corresponding to the membrane mode I fracture mode can be written as K r 0 1 v 2 T u Th sz 085 yy sin 5 preos Ha 2 o 3 v n0 n n n aan n8 gt T cos s cos 5 2 5 y Joos 2 403 SERT 2 Ky r 6 2 2 0 T v O 5 g vrsin0 v n n on n i i n 2a 2 E sin OS Bin 5 2 8 s 1 sin 2 40b where p is the shear modulus v is Poisson s ratio and E is Young s modulus The conventional finite element displacement formulation is not able to describe the first part of the displacement field corresponding to the stress intensity factor In our method this singular part of the displacement field is taken care of via the mode enrichment technique In that case the finite element formulation has to approximate the displacements for the higher order terms starting with the elastic T term The displacement field for the elastic T term varies in a linear way in x and y direction and therefore can be described exactly by the finite element displacement field However the conventional finite element formulation is not able to describe the analytical displacement fields for the other higher order terms exactly This poses a difficulty because th
23. and the analytical displacement solutions so Se JSG CyGy dV 2 34 12 I Noony UKE Kanal This also implies that matrix Sy will contain a considerable number of nonzero values because a number of elements around the crack tip are involved in the enrichment process Entering this sparse matrix together with the conventional global stiffness matrix in an equation solver would be inconvenient because it would destroy the bandwidth Therefore the corrector equations 2 33 will be solved with a Block Gauss elimination scheme as follows First equation 2 33 is written as S Aa S Ab r 2 35a ScAa SpAb r 2 35b where S Su Sp S12 F _ Ak Aa Aq Ab N page 18 B2000 User Manual 4 14 94 ra f ry fal 2 36 Working out equation 2 35a yields Aa Sr S SpAb 2 37 and substitution of equation 2 37 into equation 2 35b will result in ScSA r ScS Sg b Sp b r 2 38 or Sp SeS Sp Ab r S Sir 2 39 From this equation Ab can be calculated and Aa is obtained from equation 2 37 by substitution of Ab This technique for solving the corrector equations is already applied in the B2000 routine b2cont when dealing with the additional path equation because otherwise this path equation would completely destroy the bandwidth of the matrix Therefore the solution process discussed in this section actually implies an extension of the method that is already implemented in b2cont 2 1 1 6 A posteriori
24. are as follows step 1 step 2 step 3 Attention The first variations FVAR GLOB 0 for the mapped displacements and rotations are calculated These displacements and rotations are situated in the global displacement vector DISP GLOB O Furthermore the original and rotated orthogonal triads in the element nodes are loaded from database file TRIO iel and TRIA iel respectively These database files DISPGLOB O TRIO iel and TRIA iel are all produced by the displacement mapping routine Chapter 4 The new external load vector is calculated as F A F FVAR GLOB 0 and written to database file FRCA GLOB The final equilibrium state is calculated now by removing this external force via a load relaxation technique This option has to be used in combination with the option LMAX in order to exactly remove the residual force vector FRCA GLOB Furthermore the B2000 input file for b2contcrack has to be extended as follows LCB 2 load case 2 added PAS 1 0 Initial load value for the external load FRCA GLOB see above DPAS dA initial step for load case 1 used in the load relaxation procedure 1 0 lt DPAS lt 0 4 14 94 B2000 User Manual PAMAX 0 PBS A DPBS 0 PBMAX Xo FORCES 2 3 Additional in and output files 2 3 1 Input file mode _enrich page 29 the external forces in vector FRCA GLOB have to be released final load value for load case 1 for the nonlinear analysis
25. belonging to the finite element distribution and a set k of stress intensity factors that are coupled with the local singular fields 2 2a b and 2 3a b In this way the stress intensity factors appear as an additional set of computational freedoms In other words this technique produces the stress intensity factors as an integrated part of the set of computational freedoms The displacement functions corresponding the membrane stress intensity factors do not depend on the kind of plate theory that one is using However for the bending and tearing stress intensity factors there is a difference between the displacement solutions according to the Kirchhoff 27 and the Reissner plate theory 9 10 Because we are using a plate element that is based on a shear deformable plate theory we have to apply the analytical displacement solutions based on the Reissner plate theory For further details the interested reader is referred to 14 15 The extension of the range of applications of the method to geometrically nonlinear problems requires some careful consideration with regard to the formulation of the superposition process In the following sections we will discuss three major aspects of this approach that need clarification The first concerns the requirements that we have to impose on the size of the enriched element domain and on the continuity of the modes u and u combined Another issue which needs attention is the implementation of the local cr
26. calculates the crack tip parameters Kj Ky Ky Kpy Kgy and T with the so called mode enrichment technique This mode enrichment technique is implemented in the nonlinear Bathe element and thus enables the calculation of the crack tip parameters for linear as well as geometrically nonlinear crack problems In chapter 3 the routine for calculating a global and local error over the total finite element mesh is discussed Till now this routine can only be applied for the Bathe elements and it produces information that has to be available in order to perform adaptive mesh refinement The actual refinement of a finite element mesh is not implemented in the B2000 finite element package but is written as a stand alone program that can also be used in connection with other finite element packages An input description for this mesh generator can be found in 16 Chapter 4 deals with a routine for mapping the solutions data from one finite element mesh onto another This solution data mapping technique is very useful for geometrically nonlinear problems where the nonlinear calculation after mesh refinement can be continued without having to reload the model again from the beginning After regaining equilibrium again for the _ new grid layout the nonlinear calculation can be proceeded at the present load level This page 4 B2000 User Manual 4 14 94 technique proves to be very important for crack growth simulations where one only has to open the new crack incr
27. ction factors are needed by the QUADMESH finite element mesh generator 16 for updating a finite element mesh For this purpose output file ERRORS of routine b2error has to be copied to input file reduction for QUADMESH This format of output file ERRORS is as follows ee ERROR CALCULAT I ON ELEMENT STRAIN ENERGY ERROR 1 2 0 2 2 2 22222 2 22322 TE 2 3 eons wae ee eee ee ee ee ww x TOTAL ERROR TOTAL ENERGY 22222 PERCENTAGE ERROR REQUIRED ELEMENT REDUCTION FACTORS TO FULFILL A PERC ERROR OF 222 ELEMENT RED FACTOR 1 PEEL aes 2 ILER 3 22222 4 14 94 B2000 User Manual page 45 4 Solution data mapping 4 1 Theoretical background 4 1 1 Node projection The solution data mapping routine starts by locating the position of the new nodes with respect to the old element mesh This position is determined via a projection procedure and the solution data are mapped from the old element mesh onto the new one via the obtained finite element shape functions In the projection process the nodes of the new mesh are projected by dealing with the new elements in a sequential manner The advantage of this approach is that one can always check whether the projection points of the nodes of a new element are situated in and around one and the same element of the old mesh This way of working will prove to be very helpful in the case of a structure with two or more boundaries that are almost coinc
28. e interpolation functions of the finite element mesh try to approximate the lower and higher order terms at the same time Consequently the linear part of the finite element displacement field which should correspond to the elastic T term is affected by the higher order terms thus impairing the accuracy of the calculated elastic T term How this difficulty can be circumvented to some extend will be discussed in the next section 4 14 94 82000 User Manual page 21 2 1 2 2 Solution strategy According to equation 2 40a the analytical displacement solution u ahead of the crack tip 8 0 contains besides the elastic T term a contribution not only of the stress intensity factor Ky but also of all the other higher order terms The linear part of the finite element displacement field in this region will be disturbed significantly by these higher order terms Consequently it is difficult to compute the elastic T term from the displacement field ahead of the crack tip with a reasonable accuracy l On the other hand behind the crack tip 180 the analytical displacement solution only contains contributions from the elastic T term and from the higher order terms with even index n Hy Hs Hg etc Further there is no contribution of the displacement solution corresponding to Ky This implies that behind the crack tip the disturbance of the linear displacement field is reduced to a minimum and that the elastic T term can be computed rather accurate
29. ements obtained from each crack growth step via a load relaxation method In chapter 5 finally some special adaptations for the load processor will be discussed briefly These changes concern the fact that we want to deal with cylindrical shells and therefore have to be able to work in a cylindrical coordinate system For this purpose the external loads which are by default defined in a cartesian coordinate system have to be transformed to a cylindrical coordinate system This report does not contain any sample calculations concerning the results of the application of the developed B2000 routines The interested reader is referred to 14 15 in order to see what the capabilities of this program are 4 14 94 B2000 User Manual _ page 5 2 Crack analysis 2 1 Theoretical background 2 1 1 Mode enrichment 2 1 1 1 Basic principle For a shell or a plate the general form of the linear elasticity solution for the stresses in the immediate surrounding of a sharply edged crack does not depend on the loads and boundary conditions applied to the structure at a finite distance of the crack tip Only the stress intensity factors which are characteristic parameters in this solution depend on these conditions 5 20 It is this property that is exploited by the singular mode superposition technique The general solution for a crack in a plate or a shell in terms of the displacements u u v and stresses G G Oy Oxy tangential to the midsurface ca
30. ent identification number NFIG times 2 I5 page 36 B2000 User Manual 4 14 94 NUMNOD 7 LSQ1 3 IPOL 2 LSQ2 6 card 23 NUMNOD LS 1 LSQ2 IPOL card 24 NODEUP 1 NODEUP 2 NODEUP 3 NODEUP 4 NODEUP S NODEUP 6 NODEUP 7 card 25 IELEUP 1 IELEUP 1 IELEUP 2 IELEUP 3 IELEUP 4 IELEUP 5 IELEUP 6 card 26 NODELW 1 NODELW 2 NODELW 3 NODELW 4 NODELW 5 NODELW 6 NODELW 7 card 27 IELELW 1 IELELW 1 IELELW 2 IELELW 3 IELELW 4 IELELW S IELELW 6 IELEUP 2 IELEUP 1 IELEUP 3 IELEUP 6 IELEUP 5 TELEUP 4 NODEUP 1 NODELW 1 NODEUP 3 NODELW 3 NODEUP 2 NODELW 2 NODEUP 4 NODELW 4 NODEUP 7 NODEUP 6 NODEUP 5 NODELW 7 NODELW 6 NODELW 5 faces TELELW 3 IELELW 2 ITELELW 1 TELELW 4 IELELW 6 TELELW 5 crack tip range over which a least square fit is performed Figure 13 Definition of the input data needed for calculating the elastic T term IFIG IELE configuration number that element IELE TELE belongs to attention The element numbering may not cantian any holes and must be available in ascending order 4 14 94 B2000 User Manual page 37 3 Error calculation 3 1 Theoretical background 3 1 1 The energy norm The Z2 error estimator or energy norm is defined as 29 l 1 2 lel feg D e av i 3 1 V where o 9 9 f 3 2 In equation 3 1 D is the inverted material stiffne
31. equilibrium state is small enough this new equilibrium state can be obtained within a few iterations If this difference is too large a load relaxation procedure has to be used whereby the residual force vector resulting from the mapped solution is used as an external load 21 4 2 Interactive input commands The processor b2intp maps the solution data from an old finite element mesh onto a new one Till now the solution data mapping technique is only implemented for the 4 8 and 9 node Bathe elements The user has to enter some information concerning the names of the databases involved the last load cycle for the nonlinear calculation and the number of structural configurations This information will be asked for automatically by the program Question Enter the archival database for the new nonlinear calculation Answer The name of the archival database corresponding to the new finite element mesh has to be entered 4 14 94 Question 2 Answer 2 Question 3 Answer 3 Question 4 Answer 4 B2000 User Manual page 49 Enter the computational database from which to read the solution data The name of the computational database that was used in the nonlinear calculation with b2cont or b2contcrack has to be entered Enter the last load cycle number The solution data can only be mapped from the old finite element mesh onto the new one for the last load cycle becau
32. erentiation with respect to In our formulation we use the Green Lagrange strains which can be written as 1 T T Gj z E A u j u Aj u u 2 14 Substitution of the composite field u u gives i T G Rv 5 A u uy A 0 yy Ajt u Rv Rv Ajtu Rv Rv 215 or 1 Gi u Rv Gi u s A u i Rv Rv Aj u p v R Rv 2 16 Substituting equation 2 13 into 2 16 and making use of the orthogonality of the rotation matrix R thus RTR I yields Gy m Rv G u 5 b TRv Rv b v V3 2 17 Equation 2 17 can be rewritten as G u Rv G u G u Rv 2 18 where Gi u Rv 5 b Rv Rv b v v 2 19 page 12 B2000 User Manual 4 14 94 G describes the increase of the deformation in going from B to B It follows from the foregoing that Gi is only locally defined i e in a small area Q around the crack tip Everywhere outside this area Gi is equal to zero As long as the initial curvature and the deformation Gi u remain small in Q which is a reasonable assumption the variation of b over this area will also be small Therefore we can replace the base vectors b belonging to state B in equation 2 19 by base vectors b r b cr that are attached to a point determined by r in the region Q for which the node at the crack tip presents itself as a natural candidate Let v be specified by
33. erties of the structure surrounding the crack The solutions 2 2a b and 2 3a b are generally applicable as long as the problem remains within the realm of the theory of linear elasticity It is now noted that the function u varies with Jt while the lowest order terms of u are proportional to r For the stresses a similar qualification holds The crack tip solutions O are proportional to 1 fr while the leading term in o is a constant Due to these singular 1 Jr 4 14 94 B2000 User Manual page 7 terms the crack tip solution is dominated by u and o while away from the crack the solution field is represented by the composition 2 1a b At the crack tip r 0 all derivatives of u with respect to r are singular while away from the crack tip the solution is smooth with derivatives that are continuous in the entire remaining domain of the structure As noted 2 8 11 12 19 25 26 28 the separation of the general solution into two parts i e into a local singular field plus a regular global field can be carried over into a finite element solution scheme in an almost natural way The first step in this approach is to construct a representation of the global field u in terms of a regular finite element distribution The second step is to add to this field locally the functions u so that the completion of the solution is obtained The result of such a discretization is that the solution is represented in terms of a nodal displacement set q
34. file mode_enrich 4 3 2 Input file config 4 3 3 Input file configo Modifications for the load processor 5 1 Theoretical background 5 2 Interactive input commands 5 3 Additional in and output files 5 3 1 Input file mod _enrich 5 3 2 Input file config Examples Literature 4 14 94 44 45 45 45 47 48 49 49 50 50 51 51 51 51 51 52 53 55 4 14 94 B2000 User Manual page 3 1 Introduction The modular finite element package B2000 17 18 is extended for the analysis of cracks in a geometrically nonlinear environment For this purpose a number of routines that are already available in B2000 had to be adapted and some new numerical tools had to be developed These crack analysis tools are implemented for the 4 8 and 9 node Bathe elements 1 f This report not only serves as a user manual for the special purpose routines that were developed but also contains information about the theoretical background of these routines Chapters 2 to 5 deal with the separate routines that were adapted or developed and these chapters are written according to the following set up section 1 Theoretical background section 2 Interactive input commands section 3 Additional in and output files In chapter 2 the routine for calculating the stress intensity factors and the elastic T term will be dealt with This routine actually represents an extended version of the continuation processor b2cont and
35. g nodes of the new element is carried out with the understanding that the projection process is always started from the old element that contains the projection of the first new internal node The reason for this can best be explained by means of an example l Suppose that one has two finite element meshes for a cracked plate consisting of 4 node elements Figure 17 P the projection of the new node np which is part of the lower crack 4 14 94 B2000 User Manual page 47 old mesh sasana new element crack faces n P a es F ay lower crack A a Finite element mesh along a curved crack b Detail Figure 17 Node projection for a model with a curved crack face is geometrically situated in element e3 at the upper crack face of the old mesh Figure 17a As a consequence the mapped displacements and orthogonal triads for point P on the lower crack face will be dictated by the corresponding solution data of some old nodes at the upper crack face and it is clear that this will yield incorrect results Starting the projection process for node n from element ey the element that contains the projection of the first internal node n one will end up with element e to be the best choice for projecting node n on Due to the curvature of the crack faces however the projection point P is not situated inside element e2 Therefore P is shifted to P by setting the value of and or that is are outside the interval
36. he circumferential direction of the model and this facilitates the application of boundary conditions 2 Ifone has to deal with more than 1 structural member for instance a shell with rings and stringers one has to define these different structural members as separate element configurations For these elements the orthogonal triads which are determined as described under point 1 have to be modified in order to align them with the element surface This technique works as follows First the nodal triads are defined as in point 1 However the rotation vectors calculated in this way are not situated in the plane of the considered elements while the radial direction is not perpendicular to these elements For this purpose the orthogonal triads calculated as in point 1 have to be transformed see Figure 9 This is done according to the following rule ITRTRI IFIG 1 2 3 1 V1 2 n new for conf IFIG vi ola for conf 1 2 V12 n new for conf IFIG v otd for conf 1 3 V1 p n new for conf TFIG Vayoig for conf 1 1 V1 2 n new for conf IFIG v gia for conf 1 2 V1 2 n new for conf IFIG vo o1q for conf 1 3 v1 2n new for conf IFIG vr olg for conf 1 y original triad modified triad triad modification data for an element of configuration 2 ITRTRI 2 1 2 ITRTRI 2 2 3 ITRTRI 2 3 1 Figure 9 Transformation of the orthogonal triads when using different structural components page 34 B2000 User Manual
37. ident as for example in the case of a cracked plate First a so called internal node i e a node that is not coincident with the boundary of the model is located for each new element This procedure is followed to ensure that the projection of the initial finite element mesh usuta new element finite element boundary ee e gt physical boundary Figure 14 The projection of node ny of new element ey is situated outside the old element mesh page 46 B2000 User Manual 4 14 94 first node of each new element does not fall outside the old element mesh Figure 14 Subsequently one has to localize the node ng of the old mesh which is closest to this internal node n and also to identify an element e of the old mesh that node n belongs to Figure 15 Node n is projected onto the surface described by element e and the and coordinate of this projection point P are checked to see whether this point is situated inside element e If for instance it turns out that lt 1 0 the projection process must be repeated with the element in the old mesh that has side 1 4 in common with the previous element e see Figure 16 This technique of projecting node n and checking the size of the parameters amp and gt must be continued until and lie between 1 0 and 1 0 Figure 15 Figure 16 Values of the parameters amp and gt in and around element e The positioning of the remainin
38. ions are implemented in the B2000 processor b2error and can be selected interactively by the user After the calculation of the global percentage error the user Can enter a required percentage error that is applied to calculate the reduction factors for all the elements separately These reduction factors are used by the automatic mesh generator to build a new updated finite element model 3 1 2 Modified version of the energy norm A potential disadvantage of the error estimator described in the previous section is that it is difficult to apply to complex structures involving discontinuities These discontinuities can exist of shell intersections discrete stiffeners material or thickness jumps and or concentrated forces For such problems the model has to be partitioned at such discontinuity boundaries because of the stress and strain discontinuities that exist at these boundaries Another way to deal with problems involving physical discontinuities is to use an error estimator based on the strain energy rather than stress or strain The strain energy namely has to be continuous over these discontinuous boundaries This alternative also has the advantage of being relatively inexpensive because only one scalar quantity the strain energy density has to be dealt with instead of a tensor quantity stress or strain The continuous stress field o and strain field are obtained from the nodal stresses and strains via the shape functions Ni equatio
39. l in and output files 4 3 1 Input file mode_enrich This input file which is described already in section 2 3 1 is used only to check whether the old solution data are expressed in a cartesian or a cylindrical coordinate system i e card 1 2 page 50 B2000 User Manual 4 14 94 and 3 The mapped displacements will be defined with respect to this same coordinate system 4 3 2 Input file config This file contains the configuration data for the new finite element mesh as explained in section 2 3 2 4 3 3 Input file configo This file contains the configuration data for the old finite element mesh as explained in section 2 3 2 4 14 94 82000 User Manual page 51 5 Modifications for the load processor 5 1 Theoretical background The load processor is modified in such a way that it is possible to express the calculated nodal forces with respect to a cylindrical coordinate system This extension has to be available in order to perform a nonlinear calculation with b2cont or b2contcrack where the displacements and forces are expressed in a cylindrical coordinate system Further one has to be able to apply a pressure load to a selected number of elements For a cylindrical shell for instance the external load has to be applied to the shell elements only and not to the stringer and frame elements This is accomplished by making use of the configuration data in that the external load defined via the load proces
40. l page 23 Because the elastic T term represents a constant stress parallel to the crack faces we will only use the constant stress part to calculate this elastic T term Therefore we are interested only in the linear part of the finite element displacement field along these crack faces thus Win 5 uy u 2 42 Differentiating with respect to x yields duu 98 19 ein 3 a Ft u u J ux uy 2 43 For geometrically nonlinear problems one has to note that the crack tip region will rotate in space during deformation Therefore the in plane displacement field after rotation of the crack tip has to be considered for the calculation of the elastic T term For the strain tangential to the crack face this yields Figure 7 du u if fodu v w in OR stala S gt J J ea or 1 1 Elin 7 27 A Cu wy vp v 2 wx wi 2 2 45 The element stresses corresponding to the elastic T term can be obtained by simply multiplying these strains with the elasticity constants see equation 2 40a Now in order to establish a least square fit over these element stresses it is first necessary to choose a continuous function for the stress field behind the crack tip As already mentioned in the previous section this stress field is governed by the elastic T term and the higher order terms with even index n It follows from the previous remarks that this continuous stress function should look as follows o
41. lated in terms of small displacement gradients Actually this presents a contradiction in terms because the stresses and strains will go to 4 14 94 B2000 User Manual page 13 infinity when closing in to the crack tip This contradiction is generally accepted as a modeling error that occurs only in a very small neighborhood of the crack tip However to be consistent with the singular stress and strain field for the linearized theory the additional Green Lagrange strain Gi also has to be truncated In order to linearize the in plane strains for the rotated orientation of the crack tip zone the rotation matrix has to be determined first This rotation matrix can be calculated by using the orthogonal triads in the crack tip node for the undeformed and for the rotated crack tip zone By applying this rotation matrix to the analytical crack tip strains the in plane strains corresponding to the rotated crack tip zone can be determined From these modified strains the quadratic component in equation 2 22 can be calculated and subtracted from the total strains thus correcting the in plane strains as discussed above The truncated Green Lagrange strain Gi that has to be added to the finite element strains therefore is Gij up u3 5 vi vii 2 23 Nie The strain increments caused by insertion of the classical solutions for the crack tip behavior in the global representation of the displacements are thus obtained by a direct evaluation
42. ly along the crack faces which are situated at 180 This accuracy can even be increased by the use of quadratic 8 or 9 node elements which are able to describe the analytical displacement fields corresponding to T and Hy exactly One of the problems is that the nodal displacements along the crack face which will be used for calculating the elastic T term also have to describe the displacement field within the element they are part of Such an element covers an area for which 0 lt 180 Consequently the mode I stress intensity factor and all the higher order terms are still able to affect this linear displacement field along the crack face to some extend It is obvious that the major part of this disturbance is due to the displacement field corresponding to Ky because this field is dominant around the crack tip In our approach this problem is solved by enriching the crack tip elements However there still is a contribution of the higher order terms with even index n Hy Hg Hg etc to the total displacement field along the crack faces This means that for a general case the stress along the crack face is not constant Therefore it was decided to use the following technique for calculating the elastic T term enriched elements crack tip _ Figure 5 Refined finite element mesh around the crack tip page 22 B2000 User Manual 4 14 94 First of all the finite element mesh around the crack tip is refined by introducing a number of
43. metrically nonlinear as it is in our case the crack front may undergo a finite displacement and rotation 4 14 94 B2000 User Manual page 9 under the given load In that case the set of shape functions u should be defined with respect to the rotated geometry Consequently in the nonlinear case it is not correct to simply add the analytical crack tip solution to the global solution This can only be done after these special solutions are modified appropriately The following development is included to make this issue Clear deformed state B Figure 3 Definition of the undeformed the intermediate and the final deformed state Suppose that the shell mid surface in the undeformed state Bo is described by the mapping Figure 3 X X y amp 0 2 4 After deformation in state B the shell mid surface is described by x x y amp 2 5 where is the parameter that represents the intensity of the deformation or the time whichever is appropriate The inclusion of the evolution parameter in this discussion is purely for the sake of clarity of exposition The difference between the two states Bo and B is the displacement field page 10 B2000 User Manual 4 14 94 u u x X 2 6 The displacement u at the crack tip comprises the global field u furnished by the finite element discretization u U EO NT 8 a Q 2 7 and the singular crack tip solution u U 6 0 k G
44. mputational database that is used to save the intermediate information during the error analysis has to be entered This database is created automatically by the program and thus may not yet exist Equal distribution of the error over all the elements ENTER 1 the total volume ENTER 2 The type of error distribution has to be selected by just typing 1 or 2 see also section 3 2 1 equation 3 18 and 3 21 After entering the above information the program will calculate the local and global errors and the total percentage error This information will be displayed on the screen followed by the question to enter the wished percentage error thus TOTAL ERROR A 2 ee xxx TOTAL ENERG Y amp PERCENTAGE E R R O R Question 6 Answer Enter the required percentage error The required global percentage error has to be entered This information is used to calculate the element reduction factors that are needed by the mesh generator to refine or coarsen a finite element mesh sau a dhe A i EE LB LIB as page 44 B2000 User Manual 4 14 94 3 3 Additional in and output files 3 3 1 Output file ERRORS The processor b2error creates an output file ERRORS in which the results of the error calculation are printed This file contains for instance the strain energy per element the element errors the total energy error and percentage error and the element reduction factors These redu
45. n 3 4 Because the strain energy density U is proportional to the product of o and it is obviously inconsistent to use the same shape functions for U as for and To correct this Stanley et al 24 suggested to interpolate Ju instead of U thus Jes N NN 3 22 and this approach seems reasonable The energy norm defined in the previous section is 1 2 lel e g D o g dV 3 23 4 Continuing as above one can substitute JU for yielding Kelas selasa La BIB L SA A T emo kaksassas hka HL nal br page 42 B2000 User Manual 4 14 94 1 2 les Cf Us Aad C U fg av v 1 2 f 0 JU av 3 24 v It is obvious that in equation 3 24 the stiffness matrix D is incorporated in the strain energy densities Stanley et al 24 proved this error estimator eq 3 24 to converge to the exact error for a decreasing element size in Well n gt olea 3 25 exl Further the error estimator of equation 3 24 is a lower bound for the Z error estimator eq 3 23 It is this modified version of the energy norm that we implemented in the B2000 processor b2error It determines the energy norm after making a least square fit over the strain energy densities in the integration points These strain energy densities in the integration points are calculated in b2cont crack and are written to database file ESED ibr 0 ityp 0 3 2 Interactive input commands The p
46. n be decomposed in two Xx parts namely u u 1 6 z u r 9 z 2 1a y Figure 1 Straight crack in an infinite plate page 6 B2000 User Manual 4 14 94 o o r 8 z S r 8 z 2 1b where r and are the components of a polar coordinate system at the tip of the crack as is depicted in Figure 1 and z is the distance measured along the normal n to the midsurface of the plate The subscript 1 refers to the regular part of the solution and the subscript 2 refers to the singular crack tip solution that is dominant in a close neighborhood of the crack tip When using a plate theory that includes transverse shear deformation there is also a tearing mode that is governed by an out of plane displacement field w and by transverse shear stresses o Oz o2 For a shell in a general state of deformation the in plane crack tip solutions can be cast in the form 2 u r 0 Jr K 0 zk X 0 2 2a i 1 21 o 1 6 gt K 0 zk O 6 2 2b i VT where K i 1 2 Kp Ky k i 1 2 ky k and the particular functions X Y and depend only on the specific shell theory that is used If the transverse shear deformation is taken into account one also has w Jt Ki E 0 2 3a Gop 8 Ky 0 2 3b dr Yor r 8 The stress intensity factors in these expressions are determined by the load intensity the boundary conditions the geometry and the physical prop
47. ned in the element processor the same holds for vector 2 V1FIX fixed direction for the first in plane 3 E12 5 rotation vector only used if IFIXV1 1 and ICRSYS 2 V2FIX fixed direction for the second in plane 3 E12 5 rotation vector only used if IFIXV2 1 and ICRSYS 2 NFIG total number of element configurations I5 only used if ICRSYS 2 ITRTRI NFIG 1 1 3 transformation vector NFIG 1 times denoting how the orthogonal triads determined 3 15 in the cylindrical coordinate system have to be transformed in order to obtain the correct orthogonal triads for structural members for which the normal vectors are not aligned with the radial direction only used if ICRSYS 2 ICRCKT identification number of the crack tip 15 node NK1 NK2 NK3 NB1 NB2 NT 6 15 NK1 0 Ky is not calculated NK1 1 Ky is calculated the same holds for NK2 Ky NK3 Ky NB1 Kgy NB2 Kpyp and NT elastic T term Attention If NB1 1 and or NB2 1 and or NK3 1 and or NT 1 then also NK1 1 and NK2 1 only used if MENRI 1 CRADV 1 3 in plane vector at the crack tip node 3 E12 5 tangential to the crack faces pointing in the direction of the crack 3 only used if MENRI 1 CRADP 1 3 in plane vector at the crack tip node 3 E12 5 perpendicular to vector CRADV 1 3 3 only used if MENRI 1 NCRITM INTM number NCRITM of enriched 2 I5 elements used for the calculation of the membrane 4 14
48. new element rings inside the original crack tip elements Figure 5 Next the mode enrichment technique is applied to all the elements in this crack tip zone where the outer ring of elements serves as transition zone between the enriched and the conventional elements The enriched and transition elements are integrated with a higher order 10x10 Gauss integration scheme see section 2 1 1 7 while the conventional elements are integrated with a low order 2x2 or 3x3 integration scheme After completing the finite element calculation the element stresses corresponding to the linear finite element displacement field tangential to the crack face are determined along the element sides that coincide with this crack face The final step is then to determine the elastic T term from these element stresses This will be accomplished by making a least square fit over a number of element sides behind the crack tip The details of this technique will be described in the next section 2 1 2 3 Numerical procedure The finite element displacement solution for an element side of a quadratic quadrilateral element that coincides with the crack face can be written as Figure 6 u 34 1 8 u 1 8 us 58 1 6 241 where is one of the two parameters that are used to define the isoparametric quadrilateral element 1 1 l crack tip Figure 6 Part of the refined finite element mesh behind the crack tip 4 14 94 B2000 User Manua
49. obtained by solving this matrix equation The question remains how to choose the number of terms that have to be considered in the stress function behind the crack tip equation 2 46 in order to render the elastic T term with a reasonable accuracy Moreover it must be determined over r which elements behind the crack tip this least square fit must be established The answer to these two questions is given in 15 Summarizing one can say that for most cases one can do with a polynomial order of 2 for equation 2 46 thus G x Cyt Cix C x2 2 49 Furthermore the crack tip elements immediately surrounding the crack tip as well as the transition elements along the crack faces were shown to have a negative influence on the accuracy of the calculated elastic T term Therefore the corresponding element sides are skipped from the least square procedure In b2contcrack this can be accomplished by defining the interval LSQ1 LSQ2 of nodes over which the least square procedure is performed see section 2 3 1 input card 23 and Figure 13 2 2 Interactive input commands The new processor b2contcrack is an extension of the existing processor b2cont 4 and therefore can deal with all the interactive commands that are defined for this processor However for our purpose some special commands have been added concerning the way in which mode enrichment is applied the final load value the number of iterations for each load step the use of the pa
50. odes bordering I5 the enriched element domain used for the calculation of the bending and tearing stress intensity factors 4 only used if MENRI 1 and NB1 1 or NB2 1 or NK3 1 IPBRDBT 1 NPBRDBT identification numbers i NPBRDBT I5 of the nodes bordering the enriched element domain that is used for the calculation of the bending and tearing stress intensity factors 4 only used page 32 B2000 User Manual 4 14 94 if MENRI 1 and NB1 1 or NB2 1 or NK3 1 21 TEROT 1 2 identification number of the 2 elements 215 behind the crack tip that are used for linearizing the in plane strains gt only used if MENRI 1 22 IPROT 1 3 identification number of the3 nodes 215 behind the crack tip that are used for linearizing the in plane strains gt only used if MENRI 1 23 NUMNOD number of nodes behind the crack 4 15 tip along the crack faces used to calculate the elastic T term 6 LSQ1 LSQ2 interval of nodes behind the crack tip used in the least square procedure 6 IPOL polynomial order of the stress function used in the least square procedure 6 used only if MENRI 1 and NT 1 f 24 NODEUP 1 NUMNOD identification numbers NUMNOD I5 of the nodes along the upper crack face that are used for the calculation of the elastic T term 25 TELEUP 1 NUMNOD identification numbers NUMNOD I5 of the elements along the upper crack face that are used for the calculation of
51. of the Green Lagrange functional expressions in terms of the local components of this displacement field It is further noted that the two fields u and u are not coupled in the strain measures but only in the specific elastic energy E 0 5 Gij Cijk Gxi at least to an acceptable degree of approximation 2 1 1 4 The resulting nonlinear equations The nonlinear equations that result from the foregoing definitions can be written as F q k A 0 2 24 where represents the intensity of the loading k is the vector of stress intensity factors and q is the vector of total nodal displacements These total nodal displacements for the enriched elements are obtained by summation of the analytical displacement solution corresponding to the stress intensity factors and the underlying conventional finite element displacement field Figure 4 Thus 2 25a page 14 B2000 User Manual 4 14 94 or q q kRv 2 25b where are the nodal values of the analytical displacement solution corresponding to the stress t intensity factors are the nodal displacements determined by the underlying finite element displacement 1S formulation are the total nodal displacements _ represents the rotation of the crack tip with respect to the undeformed state and o is the vector of stress intensity factors The total displacement field can be written now as iz N amp q kRv NT E g kR v N 6 y 2 26
52. omain where T is the maximum permissible percentage error Equation 3 17 can be translated into placing a limit on the error in each element One can require for instance that the total error is equally distributed over all the elements This yields for the permissible error for each element mun _ 3 18 lel s Soom where m denotes the total number of elements According to equation 3 18 the ratio oled 3 19 em defines the elements to be refined and its value can decide the degree of subdivision or the element size needed Assuming for instance the rate of convergence to be O hP the required element size h should be h Ji 3 20 where h is the current element size Generally p is assumed to be equal to the polynomial order of the shape functions used 29 Equation 3 19 also implies that elements with an error much smaller than the required error em can be enlarged In contrast to equation 3 18 one can also require the error to be equally distributed over the total volume of the model This will yield for the maximum permissible error per element we leel S NA m 3 21 where V denotes the volume of element e and V represents the total volume of all the elements eo AMAA LBL cans 4 14 94 B2000 User Manual f page 41 Compared to equation 3 18 the maximum permissible error is proportional now to the volume of the elements 6 Both these error distribution opt
53. p precisely at the maximum load value denoted as PAMAX This means that one does not have to interpolate between the two final equilibrium states below and above PAMAX in order to obtain the solution data for PAMAX For each load step the maximum number of iterations defined in the B2000 input file as MAXIT is applied This means that although the convergence criteria are satisfied for a certain problem the maximum number of iterations will always be completed Of course this is not the case if the solutions starts to diverge and a steps size cut or a refactoring step for the stiffness matrix has to be performed This option is inspired on the fact that the bending stress intensity factors Kpy and Kgy and the tearing stress intensity factor Kyyy converge in a much slower way than the two membrane stress intensity factors Ky and Kyr However by using the maximum number of iterations the final corrections to the solutions for Kpy Kgy and Kyy are very small and convergence will be obtained 4 14 94 ILOADF NOTP B2000 User Manual page 27 The load value is fixed during each load step This option actually implies that the modified Newton method is used instead of the arc length technique When opening the new crack increment some problems occurred when using the arc length technique because the corrections to the load value became so large that it took a considerable number of load steps to reach the final equilibrium s
54. path parameter n It is characteristic for these type of procedures that each new step along the loading path is carried out by a predictor corrector scheme 21 The addition of the modes corresponding to the stress intensity factors to the regular discretization scheme however has an impact on the solution algorithm that needs some clarification The changes can best be described at the level of the corrector equations The solution of the corrector equations in each new step of the step by step solution algorithm furnishes sequentially a set of corrections Ati i 2 3 to the predicted solution 1 of the page 16 B2000 User Manual 4 14 94 equation set 2 24 This system of equations is usually a linearization of 2 24 and for the sake of simplicity we will assume that this linearization corresponds to Newton s method The governing equations on which the operations are carried out are first written as H t n 0 2 30 where F t t k eR R 1 H 2 31 oe ba The extra equation that is introduced here is a device to introduce the path parameter T that measures the distance between the solution points along the loading path Newton s method is based on the linearization H t n At H t n 2 32 which we call the Newton form of equation 2 24 The configuration t is here some approximation to the actual solution while At stands for the correction to this approximation at the prescribed value of that determines a par
55. pp 13 24 E CARDEW M R GOLDTHORPE I C HOWARD and A P KFOURI On the elastic T term in Fundamentals of deformation and fracture Proceedings of the Eshelby Memorial Symposium edited by B A Bilby K J Miller and J R Willis Cambridge University Press 1985 pp 465 476 O FOSCHI and J D BARRETT Stress intensity factors in anisotropic plates using singular isoparametric elements Int J for Numer Meth in Engrng Vol 10 1976 pp 1281 1287 J HARTRANFT Improved approximate theories of the bending and extension of flat plates in Mechanics of fracture 3 Plates and shells with cracks by G C Sih Noordhoff International Publishing ISBN 90 286 0146 5 1977 pp 45 83 J HARTRANFT and G C SIH Effect of plate thickness on the bending stress distribution around through cracks J of Math and Phys Vol 47 1968 pp 276 291 HEPLER and J S HANSEN Mixed mode fracture analysis of rectilinear anisotropic plates by high order finite elements Int J for Numer Meth in Engrng Vol 7 1981 pp 445 464 Lt nest E EE ON idee IRAR OEE TOAS LETE OAN EEA EO LAFEE OE 1 AON m MERO NENE ak ASE AAEE O ELNE ENEE EA page 56 B2000 User Manual 4 14 94 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 P R HEYLIGER and R D KRIZ Stress intensity factors by enriched mixed finite elements Int J for Numer Meth in Engrng
56. procedure b2intp Chapter 5 finally discusses the simulation of crack growth in thin walled structures and the numerical simulations described in this chapter are performed with all the processors discussed in this user manual together with the automatic mesh generator 16 which is able not only to refine or to coarsen an existing finite element mesh but also to update the mesh after crack extension 4 14 94 B2000 User Manual page 55 Literature 1 K J BATHE Finite element procedures in engineering analysis Prentice Hall Inc Englewood Cliffs New Jersey 07632 ISBN 0 13 317305 4 1982 2 S E BENZLEY Representation of singularities with isoparametric finite elements Int J for Numer Meth in Engrng Vol 8 1974 pp 537 545 3 J BLAAUWENDRAAD and A W M KOK Elementenmethode voor constructeurs 4 A 5 D J C 7 G 8 R 9 R 10 R 11 G deel 1 Agon Elsevier Amsterdam Brussel ISBN 90 10 10439 7 1973 de BOER B2000 User manual of the nonlinear analysis processor b2cont NLR CR 91470 L National Aerospace Laboratory NLR The Netherlands December 1991 BROEK Elementary engineering fracture mechanics Sijthof amp Noordhoff Alphen a d Rijn The Netherlands ISBN 90 247 2580 1 1986 BUGEDA and E ONATE New adaptive techniques for structural problems in Numerical Methods in Engineering by Ch Hirsch Elsevier Science Publishers B V ISBN 0 444 89794 1
57. rocessor b2error calculates the energy norm for a finite element model consisting of Bathe elements only A further limitation is that the elements are of one type only 4 8 or 9 node The user has to enter some computational databases and some information concerning the wished type or error calculation This information has to be entered interactively by answering some questions Question Enter the archival database used in the nonlinear calculation Answer l The name of the archival database that was used in the nonlinear calculation with b2cont or b2contcrack has to be entered Question2 Enter the archival database with the updated degree of freedom numbers shut dabesaesnledBewde LIAL 4 14 94 Answer 2 Question 3 Answer 3 Question 4 Answer 4 Question 5 Answer 5 B2000 User Manual page 43 In order to perform a least square fit over the complete model one has to solve a linear system of equations where each node has 1 degree of freedom only namely the strain energy density This database has to be created by using the B2000 input file where only the first degree of freedom for each node is free Enter the computational database used in the nonlinear calculation The name of the computational database that was used in the nonlinear calculation with b2cont or b2contcrack has to be entered Enter the computational database used for the error calculation The name of the co
58. s means that the formation of the stiffness matrices for the enriched elements is rather time consuming due to the high number of integration points and the evaluation of the additional stiffness terms in these integration points A posteriori application of the mode enrichment technique is thus an effective means to reduce the volume of the computation because the higher order integration scheme has to be applied only once namely after the nonlinear deformation of the structure is completed This technique is implemented in the B2000 processor b2contcrack via an interactive command ENRI see Section 2 2 Furthermore another technique is implemented where mode enrichment is applied for the complete loading history while the calculation is repeated once for the final load step without enriching the crack tip zone via the interactive input command REMO These final results can be used to perform an error calculation because the displacements are now governed by the finite element formulation only This also has the advantage that after the solution data is mapped from an old finite element mesh onto a new one with an extended crack the old crack tip region in the new finite element mesh does not have to be enriched when relaxing the residual loads that keep the new crack increment closed 2 1 1 7 Numerical integration In finite element analysis the integrals involved in the equilibrium equations are evaluated numerically The type and
59. se the updated triads in the element nodes are available only for this last load cycle Note that the displacements have to written in branch local coordinates i e DISP 1 icycle In version 1 77 of B2000 this branch displacement vector has to be generated by running the b2ebv processor Is the number of configuration larger than 1 enter 0 for no or 1 for yes If one is dealing with a model consisting of different structural configurations b2intp has to take care of the fact that the mapped solution data for each new element are obtained from old elements of the same structural configuration This is especially important for the orthogonal triads in the element nodes because the definition of these triads depends on the shape of the separate elements that the same node is part of This information is used to map the solution data from the old finite element mesh onto the new one The mapped displacements are written automatically to the new branch displacement vector DISP 1 0 and to the global displacement vector DISP GLOB 0 This global displacement vector is used as an initial estimate to restart the nonlinear calculation Further the nodal triads in the original and the deformed state are written to TRIO iel and TRIA iel respectively These initial displacements and orthogonal triads are entered in the nonlinear calculation by using the option NMESH see section 2 2 4 33 Additiona
60. sor is always applied to elements of configuration 1 only 5 2 Interactive input commands No additional interactive input commands have to entered 5 3 Additional in and output files 5 3 1 Input file mode_enrich This input file which is described already in section 2 3 1 is used only to check whether the calculated nodal forces have to be expressed with respect to a cartesian or a cylindrical coordinate system i e card 1 2 and 3 page 52 B2000 User Manual 4 14 94 5 3 2 Input file config This file contains the configuration data for the finite element mesh as explained in section 2 3 2 The external load that is applied via the load processor is restricted to elements of configuration 1 only 4 14 94 f B2000 User Manual page 53 6 Example Examples concerning the routines discussed in chapter 2 to 5 can be found in 15 In chapter 2 of this reference some examples concerning the mode enrichment process are dealt with These problems only make use of the processor b2contcrack described in chapter 2 of this manual Next the calculation of the elastic T term is discussed and applied to several classical and geometrically nonlinear crack problems in chapter 3 of 15 In chapter 4 adaptive mesh refinement is dealt with For the sample calculations in this chapter the B2000 processor b2contcrack is applied together with the error calculator b2error and for geometrically nonlinear problems the displacement mapping
61. ss matrix and V denotes the volume of the finite element model The finite element stresses are represented by g equation 3 2 while is the still unknown close approximation of the exact continuous stress field One way to approximate the exact stress field is by assuming the same interpolation functions for this stress field o as for the finite element displacement field u This implies that the interpolation polynomials used for the approximate stress field o are one order higher than the ones used for the finite element stresses o which are described by the derivatives of the interpolation functions for u The nodal values of the stress field o are determined via a least square fit with respect to the finite element stress field o So if u Nu 3 3 we now assume that co No l 3 4 page 38 B2000 User Manual 4 14 94 where o represents the discrete nodal values of the still unknown continuous stress field 29 Next an error E in the stresses is introduced as 3 E o 9 av y Working out equation 3 5 yields Substitution of equation 3 4 into equation 3 6 yields for each element P Ys f T ES J oN Nio dV J oN odv 1 f gaV 24 i 2 or 1 T T Es 59 B o 9 fe tA where T B Ni N av v T f JN gav A 5 gav V Summation of the contribution Eg over all the elements yields for the total error E 50 Bo olf A This error E has to be minimized with respect to o
62. tate However by fixing the load value it turned out that the final equilibrium state could be obtained within 1 or 2 load steps This option is often used in combination with the option ITMAX An imperfection is introduced in order to avoid the existence of bifurcation points This option is used to pass bifurcation points in a smooth way by using the continuation processor In the first load step the total applied external load is F dA F A F 1 14 1 2 2 while for the rest of the loading history the external load is denoted as Fa A d F F F 1G 1 1o 1 t4272 The subscripts without brackets denote the load case number while the subscripts within brackets represent the load step number Attention When using this option the B2000 input file for p2contcrack has to be extended as follows LCB 2 load case 2 added PBS gt load value for the imperfec tion load DPBS 0 PBS may not increase dur ing the loading history PBMAX 0 FORCES case 2 imperfection forces have to be defined as case 2 page 28 NMESH B2000 User Manual 4 14 94 The nonlinear calculation is proceeded again after mapping the solution data from the old mesh onto the new one This option implies that the equilibrium equations start from a deformed state where this state does not represent an equilibrium state for the new finite element model The succeeding steps involved when using this option
63. th following technique imperfections and information to denote whether the calculation is restarted after applying the solution data mapping routine Command Description MENRI Mode enrichment is applied in the crack tip elements and the stress intensity factor s and eventually the elastic T term are calculated for at least 1 load step page 26 ENRI REMO LMAX ITMAX B2000 User Manual 4 14 94 Mode enrichment is applied only for the final load step This way of working can be used if one is only interested in the final values for the stress intensity factor s and elastic T term It may cause a considerable reduction of the CPU time as compared to the application of mode enrichment for all load steps because a lower order integration scheme is applied throughout a substantial part of the nonlinear calculation see section 2 1 1 6 This command can only be used in combination with MENRI Mode enrichment is applied during the whole loading history However for the maximum load value an additional calculation is performed without applying the mode enrichment technique This additional information can be used to perform an error calculation chapter 3 or to simplify the opening process of the new crack increment It removes the necessity of enriching the old crack tip region in the new finite element mesh as the crack extends This command can only be used in combination with MENRI The nonlinear calculation will sto
64. ticular point on the loading path In matrix form equation 2 32 can be written as Sr Sp F Aq f T _ Siz S22 Fial Ak lhl i 2 33 n n no A 0 where Sj is the global stiffness matrix corresponding to the conventional finite element displacement formulation Si is the global sparse stiffness matrix containing the coupling terms between the conventional displacement formulation and the analytical crack tip displacement solution So is the global stiffness matrix containing terms due to the analytical crack tip displacement solutions only 4 14 94 B2000 User Manual page 17 F is the vector of nodal forces corresponding to the conventional displacement field Fk is the load vector corresponding to the analytical crack tip displacement solutions n is the tangential vector to the solution field for the current load nk is the stress intensity part of the total tangential vector Ng is a constant that defines the auxiliary surface needed in the continuation procedure Ak is the correction of the stress intensity factor s Aq is the correction of the nodal displacements AX is the correction of the load Step size fi is the vector of residual forces corresponding to the conventional displacement field and fk is the vector of residual forces corresponding to the analytical crack tip displacement solutions As already stated above Sj contains the coupling terms between the conventional displacement formulation
65. w 6y yotapues jo ssew oSI1JIO2 ds fo i tIS3uy x 1Ix ZS19 g8z 1oo3 s 3euTpiIoooS BuTueu 1 5 969zZ SySb3 b 5b3auy bpo3ur TTZ 0603 bquy oquy uu sassauyotyq 2 z 53ur zo3qurv b 060U 1 63uy 19 uy Z SPSpu watee 9 s d Jo z qumu a zeya ztud teua tyud pr edut zt edut pe tduy 3 9 F8LE T pweToyur EZ Tduy pz ez 7z rz pA A ZA TA px Ex Zx TX 3 9 380 0 1I0204I GT PT ET ZT TT 000 Z 000 A o00e x Tear ssum by saybyem 27Jpo ds 2 b Gauy poquygogquy Zoquy Toquy Jabaquy u Zabayquy yo Ttdwy Z o0 Yy e Teaz Apo Tdwy z 0 y e Tear AYpoTTdut watee Byjuy 2 A x wauw euTynozqns pus dOs s Z BVTIM 10 44 ITIM 3SQN3 s xZ 3120 O a Z O9TIM ZaLIYM uranoaqns pua poquz 2 poqut A pO3uT x g 3uy 2 33u1 g33u 1 x goquz Z zaqut A z 3u 1 x t 3u1 2 tT93ur1 A T93u1 x 3u u T penb jo s pou i u Io5o5 S 3BUIDI009 ost 37 IOTOO a g O01J 01J cZ 0TJ ET 9 O1J Z8 01J ZZ TOUS AT E OTIS TEs 9 JEOTI lZ SE OTT UL e E OUT OE E OTT 02 3 OUI 0T 3 E91 629 TT3HS 8 BOWAGE 6 00 nZ BVFIM p bauz poqut g 6quy Eoquy 2 6quy zoquy 1 63uzaToquy o00 2 000 4 000c x Te 1 p Buz posuy eoquyz zoquy T237 136377 6quy 2 A X TALINM uyqnoiqns ww ew pus SAILIING Z 3TIM Zels Z

Download Pdf Manuals

image

Related Search

Related Contents

Manual de Instalação DuoBlock Fan/Titan modelo 2009  ➊SET ➋OK + - - Leroy Merlin  Instrukcja obsługi EN  Samsung SCX-4720FN Felhasználói kézikönyv  Employer Internet User manual Local Government Pension Scheme  Bedienungsanleitung  Philips Micro SD cards FM02MD35B  USER MANUAL MANUALE D'USO MANUEL D  CONFIGURACION BH 751  Moduli di sicurezza 安全继电器 Relé de segurança Relé de  

Copyright © All rights reserved.
Failed to retrieve file