CTU Prague & TU Eindhoven Mechanical Engineering, Mechanics of Solids

Experimental, numerical procedure to establish Paris’ law for verification of user subroutine for prediction of crack propagation in Abacus Internship abroad 4W708 MT 07.06

L.A.A. Beex

Supervision: M. Španiel J. Jurenka M.G.D. Geers

Prague, February 2nd 2007

Table of Contents 1. Introduction ...................................................................................................................................3 2. User subroutine for crack propagation ..............................................................................4 3. Working procedure .....................................................................................................................5 4. Script for fast computation of J-integrals for different geometries.......................10 4.1 Geometry ..........................................................................................................................................10 4.2 Material behavior.............................................................................................................................12 4.3 Mesh .................................................................................................................................................12 4.4 Boundary conditions and loads ........................................................................................................15 5. Experiments .................................................................................................................................16 5.1 Specimens .........................................................................................................................................16 5.2 Equipment.........................................................................................................................................16 5.3 Problems with draw bench and micro camera .................................................................................17 6. Results ............................................................................................................................................19 7. Conclusion ....................................................................................................................................26 8. Symbols ..........................................................................................................................................27 9. Literature ......................................................................................................................................28 10. Acknowledgements .................................................................................................................29 11. Appendices .................................................................................................................................30 11.1 Appendix A: Geometries of specimens ...........................................................................................31 11.2 Appendix B: Theoretical results for fatigue crack propagation .....................................................33 11.3 Appendix C: 3D & 2D FE-model of specimen in Abacus...............................................................34 11.4 Appendix D: Python script for fast computation of J-integrals for different geometries................35 11.5 Appendix E: Crack tips with and without rosette mesh configuration............................................51 11.6 Appendix F: Specimens & Equipment ............................................................................................52 11.7 Appendix G: Images from micro camera........................................................................................55

2

1. Introduction In a lot of applications it is important to know how cracks may propagate through material under cycling loads and how many cycles are needed to obtain a certain crack length (see Figure 1.1). For airplane wings, for example, it’s necessary to know after how many cycles, and thus flights, the critical places where cracks may occur have to be checked. It’s also important to know how a crack propagates so that this can be kept in mind during the design stage. A good example is to arrest a crack by making a hole next to a critical point where a crack is likely to occur. A user subroutine for FE-program Abacus is created by the Czech Technical University in Prague to predict crack propagation under cycling loads in metal specimens for which plane strain conditions. It’s only meant for high cycle fatigue problems so the use of linear elastic material behavior is sufficient. The problem is that the validity of the script isn’t proven yet. Therefore the aim of this research is to test if it is correct. To predict crack propagation and the number of cycles necessary to obtain a certain crack length, next to the Young’s modulus and the Poisson’s ratio, also Paris’ law has to be known. A combination of experiments and numerical analyses are needed to define Paris’ relation. The crack propagation per number of cycles will be determined by experiments, while finite element models of the test specimens have to qualify the stress intensity factor. Once Paris’ law is known, it can be implemented in the user subroutine to determine the number of cycles matching the crack path. The user subroutine is valid if the crack path determined by the subroutine and the matching number of cycles is the same as those of the test specimens. In this report, at first the user subroutine will be shortly explained. Subsequently, the working procedure to verify the user subroutine will be discussed. Afterwards, a Python script will be discussed that is written in Complete Abacus Environment (CAE) so that with little efforts the stress intensity factor can be determined for different geometries. After that, the experiments are discussed and the results of the experiments and numerical solutions will be combined to determine Paris’ law. Finally this will be used to determine the crack path with a 2D-model combined with the subroutine and this will be compared with the real crack path in the test samples.

Figure 1.1: Crack propagation through a metal specimen in one of the tests.

3

2. User subroutine for crack propagation To predict crack propagation in 2D geometries under cycling loads, the Czech Technical University in Prague has created a user subroutine, written in Python in CAE. With the subroutine, the crack path and the number of cycles, necessary to obtain this crack path, can be determined. It can only be used for 2D geometries for which plane strain conditions are valid. Besides the geometry and the linear elastic material parameters, also the locations of crack tips, loads, boundary conditions, the wanted total crack length and number of steps, to divide the total crack length in partial crack lengths, have to be defined. Using this information the user subroutine automatically creates a finite element model, including rosette mesh configurations (see Figure A in Appendix F) at the given crack tips which are necessary to correctly calculate the stress intensity factors (KI and KII). If the stress intensity factors are calculated, they have to be averaged according to Formula 2.1. The newly obtained stress intensity factor can then be implemented together with the partial crack length in Paris’ law (See Formula 2.2) to obtain the number of cycles, necessary to develop the partial crack lengths. In this case, the load changes between 0 and a given load so that is why it’s allowed to write K instead of ΔK in Formula 2.1.

K = K I2 + K II2 da m = C (K ) dN

[2.1] [2.2]

In Paris’ law, da/dN is the crack propagation per number of cycles. This is a function of the stress intensity factor and parameters C and m. These parameters are supposed to be material independent, but in practice they are often also dependent of geometry, load and frequency. Attention should also be paid to the fact that Paris’ law is only valid on a specified domain. In Figure A in Appendix B this domain is marked with II; only here Paris’ law can be determined with a power law. The maximum tangential stress criterion is used to determine the direction of the crack path. The criterion is based on the idea that the crack propagates in the direction as in which the shear stress equals 0. In Formula 2.3 the near-crack-tip shear stress is given for a homogeneous, isotropic linear elastic material. If this equals 0, the direction of crack propagation can be given with angle θ in Formula 2.4.

τ rθ =

⎛θ ⎞ cos⎜ ⎟[K Ι sin (θ ) + K ΙΙ (3 cos(θ ) − 1)] 2πr ⎝2⎠ 1

⎛ 3K 2 + K 4 + 8 K 2 K 2 ΙΙ Ι Ι ΙΙ 2 2 ⎜ K Ι + 9 K ΙΙ ⎝

θ = cos −1 ⎜

⎞ ⎟ ⎟ ⎠

[2.3]

[2.4]

Combining the direction of the crack path and the partial crack length, the user subroutine can make a new geometry with the new partial crack path. After this, the user subroutine repeats the whole procedure as many times as necessary to reach the wanted crack length. The number of cycles necessary to obtain the wanted total crack length is simply the sum of number of cycles of all partial crack lengths.

4

3. Working procedure

To check if the user subroutine is correct, first the unknown material parameters in Paris’ law have to be established (C and m in Formula 2.2). This has to be done by a combination of experiments, with specially designed compact tension samples (see Figure A and B in Appendix A), and numerical simulations of the experiments. The crack length per number of cycles will be determined by the experiments while at the same time numerical simulations of the experiment have to define the J-integral (See Formula 3.1). In this Formula W is the strain energy, ti are stress components, and ui are deformations. ⎛ ∂u ⎞ J = ∫ ⎜⎜Wn1 − t i i ⎟⎟dΓ ∂x1 ⎠ Γ⎝

[3.1]

The J-integral is a contour integral and relates energy release with crack growth and it’s a measure for the deformation of a crack tip. In our case, the material behavior is linear elastic and than the J-integral can be related to the stress intensity factors. For plane strain conditions, the relation is as follows:

(1 −ν ) (K J= E 2

2 Ι

)

+ K II2 +

1 +ν 2 K III E

[3.2]

With the J-integral, the first stress intensity factor can be easily calculated for 2D plane strain models in mode I according to Formula 3.2 because in case of mode I loading, KII and KIII both equal zero. A Python script is written so that with little efforts the J-integral, and thus the first stress intensity factor, can be calculated for different geometries that are loaded in mode I (see Figure 3.1).

Figure 3.1: The three different modes of crack tip deformation. From left to right: mode I (tension), mode II (in-plane shear) and mode III (out-of-plane shear).

If the stress intensity factor and the crack length per number of cycles are both implemented in Paris’ law, the undefined material parameters in Paris’ law can be determined. Subsequently Paris’ law can be implemented in the user subroutine and after that, new experiments have to be done to determine if the user subroutine predicts the correct trajectory of the crack and the correct number of cycles to

5

establish the crack path. The test samples in these new experiments have almost the same geometry as in the former experiments only in this case, a hole is drilled in the samples. The new specimens are not symmetric, due to this hole, and therefore the direction of crack propagation won’t be straight. So not only the user subroutine is validated for the predictions of the correct number of cycles but also for the prediction of the correct crack path. Before the unknown parameters in Paris’ law can be defined, first has to be proven that plane strain conditions are valid for the used test samples, because the user subroutine can only be used for 2D plane strain models. This is done by modeling the test sample in Abacus in 3D, 2D with plane strain conditions and 2D with plane stress conditions (See Figure C in Appendix C). One has to assume that the 3D model of the test sample gives the correct value of the J-integral, because it can’t be determined in another way. The J-integral can also analytically be calculated for plane strain models. If this matches the J-integral determined by the 2D plane strain model in Abacus, one can assume that it’s allowed to use plane strain conditions for the test samples. In the experiment to establish Paris’ law, two pins are placed in the holes of the test specimen. A force is prescribed for one of the pins, while the other one is kept in place. The specimen is able to rotate around the pins. This means that in the finite element models of the specimens an extra node is created in the two holes and the xand y-movement of the nodes at the inside of each hole are constrained to those of the extra node by couplings. The lower extra node is fixed in x- and y-direction, while the upper extra node is fixed in x-direction and the force is prescribed in y-direction (See Figure 3.2). This means that the rotation of both holes is free. These conditions match those of the experiments, except that for the upper hole only the upper nodes should be constraint to the extra node and for the lower hole only the lower nodes should be constraint, but because this is complicated to establish and research pointed out that the value for J-integrals only differs 5%, this isn’t done in the numerical analysis. Figure 3.2: 2D FE-model of the test sample, including boundary conditions, load and constraints.

Another point of attention is the definition of the J-integral in Abacus which is calculated with virtual crack extension method. This method relates energy release rate to crack extension and is more accurate than original J-integral calculations and not so sensitive for mesh size of the rosette tip In Abacus a number of contours for calculation of the J-integral can be defined. The first contour consists of those elements directly connected to the crack tip nodes (red in Figure 3.3). The next contour consists of the ring of elements that share nodes with the elements in the first contour as well as the elements in the first contour (white in Figure 3.3). Each subsequent contour is defined by adding the next ring of elements that share nodes with the elements in the previous contour. The J-integral is supposed to be path-

6

independent, but this has to be established. This is done for all the FE-models. For the 3D model the J-integral changes over the thickness. Figure 3.4 shows how the elements over the thickness are defined for the 3D model. For a force of 10kN, Emodulus of 209 GPa. and Poisson’s ratio of 0.3, all the results can be found in Table 3.1.

Figure 3.4: Definition of the elements over the thickness for the 3D model.

Figure 3.3: 4 Different contours of the J-integral in a rosette mesh configuration in Abacus.

Element# (over thickness)

J-integral *103 [N/m]

Contour 3D model

1 1 1 1 2 2 2 2 3 3 3 3 -

1 2 3 4 1 2 3 4 1 2 3 4 2D model – Plane strain 1 2 3 4 2D model- Plane stress 1 2 3 4

1.20 1.22 1.22 1.22 1.46 1.49 1.49 1.49 1.73 1.74 1.75 1.75 0 1.56 1.56 1.56 0 1.72 1.72 1.72

Table 3.1: Results of J-integral for 3D-model, 2D-model with plane strain conditions and 2D model with plane stress conditions.

7

The numerical results point out that for all cases, the J-integral doesn’t change after the third contour. This means that for all the next calculations the J-integral, matching the third contour, can be used. Another point that can be concluded from the results is that for the 3D model the J-integral grows inwards the specimen. This can be explained after careful investigation of the crack tip. The J-integral consists of a energy term and a stress component in a direction multiplied with the derivative of the deformation in the same direction to the x-direction (see Formula 3.1). The elastic energy changes over the thickness and is bigger on the sides than in the center. This means that the energy term isn’t the cause that the J-integral is bigger in the center. If we look at du1/dx1, we see that the term itself is relatively small and that it hardly changes over the thickness. Also du2/dx2 is small and doesn’t change. On the other hand du3/dx3 is large at the sides of the specimen but in the middle of the specimen it’s really small, so this is the term responsible for the change of the J-integral over the thickness. It’s also clear that the J-integral of the 2D model with plane strain conditions is smaller than the one of the 2D model with plane stress conditions. The reason for this is that a force is prescribed for the models and not a displacement. A plane strain model is stiffer than a plane stress model because a plane stress model can move in the third direction. So if a load is prescribed, the plane strain model will not deform as much as a plane stress model and thus is the J-integral of the plane stress model larger. The 3D model is assumed to give the correct J-integral. If this is averaged the J-integral is (1.22+1.49+1.75)*103/3=1.49*103 N/m. The J-integral matching the 2D plane strain model is quite similar to this one with only 4% aberration. This is better than the plane stress model because this one has an aberration of 13%. Now the analytic result will be considered so that can be said for sure that it’s allowed to use a plane strain model for our specimens. The following formula is derived by Srawley for these tests: ⎡⎛ a ⎢⎜ 2 + W ⎢⎝ P KΙ = ⋅⎢ B W ⎢ ⎢ ⎣

2 3 4 a ⎛ a ⎞ ⎞⎟ ⎤ ⎛a⎞ ⎛a⎞ ⎞⎛⎜ − + 0 . 886 4 . 64 13 . 32 14 . 72 5 . 6 + − ⎜ ⎟ ⎟⎥ ⎜ ⎟ ⎜ ⎟ ⎟⎜ W ⎝W ⎠ ⎠⎥ ⎝W ⎠ ⎝W ⎠ ⎠⎝ ⎥ 1.5 a⎞ ⎛ ⎥ ⎜1 − ⎟ ⎥ ⎝ W⎠ ⎦

[3.3]

In this formula, P is the applied force, B is the thickness of the specimen, W is width from the point where the force is applied to the end of the specimen and a is the crack length from the point where the force is applied to the crack tip. All this is shown in Figure 3.5.

8

Figure 3.5: 2D model of the specimen with all the parameters necessary to analytically determine KI.

In our case the applied force is 10 kN, W equals 0.059 m, a equals 0.0125 m and B is 0.01 m. If we implement this in Srawley’s formula, a stress intensity factor is obtained of 17.97 MPa*√m. Subsequently, formula 3.2 is used to determine the J-integral. This equals 1.41 kN/m. So this is an aberration of 10% with the plane strain model and 20% with the plane stress model. Concluding one can say that’s allowed to use the plane strain model because the results compared to the 3D model only differ 4% and compared to the analytical results differ only 10%.

9

4. Script for fast computation of J-integrals for different geometries

A lot of numerical simulations have to be done to determine the first stress intensity factor for many 2D geometries who are loaded in mode I and for which plane strain conditions are valid. To do this in an efficient way, a Python script is written in CAE which calculates the J-integral (See Appendix D). The idea of the script is that one only has to determine a 2D geometry, where the crack tips are located, the linear elastic material parameters and the boundary conditions and loads and than the script will automatically calculate the J-integral. 4.1 Geometry First of all, it has to be clear that all the dimensions of the geometry have to be given in micrometers. This is done to prevent problems with rounding off errors that may occur in Abacus at a magnitude of 10-6. A consequence of this is that applied forces have to be given in microNewtons and the Young’s modulus in MegaPascal. This means that the evaluated stresses will be calculated in MegaPascal and the J-integral in Newton/meter, exactly what is wanted. The input of the outer geometry in the script consists of the implementation of coordinates of points. The script connects these points with straight edges (see Figure 4.1.1). This can be easily done if all the points are in one list in clockwise order. Only for the last edge it’s not trivial because at the end of a list it’s not possible to skip to the first member of the list. Lists in Python can be seen as circles except at the end one can only move from the first member to the last but not from the last to the first (see Figure 4.1.2). So an extra line has to be implemented in the script for the last edge.

Figure 4.1.2: Symbolic representation of a list in Python. One can see a list as a circle and one can skip from member to member and back, except it’s not possible to skip from the end of the list to the beginning.

Figure 4.1.1: Geometry is created by straight edges.

10

The fact that the outer geometry is subscribed by only straight edges has the disadvantage that if the geometry contains an arc, this has to be defined by a number of points. The advantage, though, is that it’s transparent for the user how geometries are created and that every possible geometry can be made. The user also has to determine which of the given points are crack tips. After the script has created all the edges, it will delete the two edges that are connected to a crack tip (see Figure 4.1.3). Subsequently, four points are created on the same lines as the deleted lines, two at a distance of 1.5mm and two at a distance of 0.15mm from the crack tip. These points define the start and end point of two circles which have to be created for the rosette mesh configuration. The crack tip is the center point of these circles. The inner elements of the rosette mesh configuration will be located at the smallest circle and the outer elements of the rosette mesh configuration are located at the biggest circle. Also two new lines have to be created from the two new points closest to the crack tip to the next two points of the geometry. The actual sequence can be seen in Figure 4.1.3, 4.1.4 and 4.1.5. First the inner circle has to created, than the two new lines and after that the outer circle can be added. If a point of the outer geometry, but not a crack tip, is at a distance of less than 1.5mm from the crack tip, the outer circle of the rosette mesh configuration will not be at a distance of 1.5mm but at a distance of the point.

Figure 4.1.3: Edges near crack tips are removed. Arcs are created around crack tips.

Figure 4.1.4: Arcs are connected by straight edges to surrounding edges

11

Figure 4.1.5: Part is created

It may be clear that only the outer geometry can have crack tips. The inner geometries can consist of an infinite number of circles and rectangles. The circles have to be defined by the center point and a point on the circle. To create a rectangle, the coordinates of the four corner points have to be defined in clockwise order. When the geometry is made and the thickness is determined, the part can be created. The script will also check the distances between crack tips. If two crack tips are within a distance of 3.1mm from each other, the script cancels its work and sends out an error message to the user. This has two reasons; at first, the J-integrals can be regarded as meaningless because the two stress fields of the crack tips can influence each other. At second, it is impossible to make a decent mesh, because the two rosette crack tips are placed over each other and the mesh generation will go wrong or the space between the rosette crack tips is too small to put a normally shaped element between them. 4.2 Material behavior The script can only be used for linear elastic materials. It is not necessary to use more complicated material behavior because Paris’ law can only be used for linear elastic materials because the plastic zone at the crack tip has to be negligible small. The material parameters that have to be defined are therefore only the Young’s modulus and the Poisson’s ratio. 4.3 Mesh All the elements in the mesh are quadrilateral quadratic. This means that the elements have four corner nodes and four nodes at the ribs, which makes the shape functions quadratic and therefore more accurate than linear shape functions. The nodes have 3 degrees of freedom; x- and y-direction and z-rotation. The mesh size has to be determined by the user.

12

The most special and important parts of the mesh are the rosette mesh configurations. These parts have the purpose to describe the stress field near the crack tip properly. The stress gradients near the crack are namely very high and therefore the mesh has to be denser towards the crack tip. The rosette mesh configuration enables the mesh to go from relatively large, outside the rosette crack tip, to very small near the crack tip (see Figure 4.3.1). This is achieved by seeding the straight edges biased that beacon the rosette crack tip: there are five rows of elements and the outer elements are three times larger than the inner elements. In Figure 4.3.1 this can be found but for four rows of elements. The rosette mesh configuration made by the script at this point has still an empty circle with a radius of 0.15mm at the crack tip. To complete the mesh at the crack tip, the inner nodes at the crack tip are placed over each other, exactly at the crack tip (see Figure 4.3.2). These nodes are constraint to each other in x- and y-direction and zrotation because they all represent the same material point. For materials with a high Young’s modulus this isn’t necessary because the value of the J-integral differs negligible. On the other hand, it’s necessary for materials with low Young’s modules, so the constraints are always in the FE-model created by the script. Finally the middle nodes of the inner elements have to be moved to a distance of a quarter of the element length (see Figure 4.3.3) to describe the stress singularity at the crack tip.

Figure 4.3.1: Mesh is created including rosette tips

Figure 4.3.2: Inner nodes are placed over each other at the crack tip and the middle nodes of the inner elements are placed at a distance of a quarter of the element length.

13

Figure 4.3.3: A collapsed inner element of the rosette mesh configuration with midnodes at a quarter of the distance of the element.

The configuration of the inner elements of the rosette mesh has to be like in Figure 4.3.3 because only in this case the derivative of the stress of the crack tip can be infinite, which is the case for cracks. This can be easily shown for a 1D element. The case that the midnode is located at a quarter of the distance of the element size (see Figure 4.3.4), is worked out.

Figure 4.3.4: 1D element; x1=0, x2=L and x3=¼L.

A choice can be that x1 and x2 are respectively 0 and L, x3 will be open for the moment. In this case, the position, the displacement (x) and the stress (t) can be described with the next formulas: 1 x = ξ (ξ + 1)L − ξ 2 − 1 x3 2 1 1 t = ξ (ξ − 1)t1 + ξ (ξ − 1)t 2 − ξ 2 − 1 t 3 2 2 [4.3.1] dt dt dt dξ = = dξ dx dx dξ dx dξ

(

)

(

)

If we now consider the special case that x3 equals ¼L, the following results are obtained: 1 1 1 4x 2 x = ξ (ξ + 1)L − ξ 2 − 1 L = (ξ + 1) L ⇒ ξ + 1 = L 2 4 4 dt dt dt dt dξ = = dξ = dξ dx xL dx dξ dx [4.3.2] dξ

(

)

dt =∞ dx ξx ==0−1

14

This is exactly wanted because the derivative of the stress should be infinite at the crack tip in linear elastic fracture mechanics. 4.4 Boundary conditions and loads Boundary conditions and loads can be prescribed in three different ways for FEmodels created with the script. In all the cases an extra node is created for which the boundary conditions and loads are prescribed. Subsequently, the extra node is constrained in x- and y-direction to other nodes of the FE-model. These can be nodes at the circle or rectangle of an inner geometry or nodes at an edge of the outer geometry (See Figure 4.4.1, 4.4.2, 4.4.3).

Figure 4.4.1, 4.4.2 and 4.4.3: the three possible ways to apply boundary conditions and loads, an extra node is created in the center of a circle, a rectangle or in the middle of an edge and boundary conditions or/and loads are applied to that node. The nodes at the circle, rectangle and edge are constrained to the extra node in x- and y-direction.

15

5. Experiments

The first experiments have to establish the crack propagation per number of cycles. A cycling load will be applied to the compact tension specimen in the experiments. The crack propagation will be recorded with a micro camera while the number of cycles is counted by the draw bench. With the crack propagation per number of cycles and the stress intensity factors from the numerical analysis, the unknown parameters in Paris’ law can be determined. Subsequently, Paris’ law has to be implemented in the user subroutine and it can be checked on the hand of experiments with a slightly different geometry. 5.1 Specimens The specimens used in the experiments to determine the unknown parameters in Paris’ law have to fulfill a number of requirements. First of all, they have to be made of metal. The steel is CSN EN 100083-1+A1. The Young’s modulus and Poisson’s ratio are respectively 209 GPa and 0.3. The second requirement of the specimens is that plane strain conditions have to be valid. Therefore the thickness is set to 10 mm. The geometry of the samples can be found in Appendix A, Figure A. Another requirement is that it has to be exactly known where the crack starts growing. To accomplish this, a notch is added to the specimen. This notch has a V shape over the thickness to make sure that the crack starts growing in the middle of the specimen and not on a side of the specimen because the crack has to propagate symmetrically through the specimen. Also two holes, in which pins can be placed, are drilled in the specimens to be able to apply a cycling load (see Appendix F, Figure D). The holes are slightly larger than the pins what makes it possible to easily place them in the holes. Another consequence is that the specimen is not fixed because it can still rotate around the pins. So it is only fixed in x- and y-direction. The last requirement is that the specimens are easy to produce, because of costs and time. Therefore the geometry is kept as simple as possible. This is also the reason that the outer geometry mainly consists out of an arc. The specimens, used to check the user subroutine after Paris’s law is implemented, have almost the same geometry as the first specimens, only an extra hole is added (see Figure B in Appendix A). This hole makes the specimens asymmetric. The crack will now propagate asymmetrically which makes it possible to check the user subroutine not only on the correctly predicted number of cycles but also on the correctly predicted crack path. 5.2 Equipment The necessary equipment to perform the experiments consists of a draw bench with control unit and a micro camera with light source and hardware for visualization. The used draw bench is a 30 year old Heckert UFP 400 with a range of 80 kN (see Figure A in Appendix F). The maximum frequency of the machine is 10 Hz. In the experiments the lower boundary of the force is always 1 kN, because if the force is set on 0 N it overshoots. This may lead to problems with the hydraulic actuation so this has to be prevented. The source of the overshoot may have several reasons. It can be

16

caused by bad sensors and a low sampling frequency. Another possible explanation is the cold weld that appears between the pins and the specimen. Due to this, the pins are stuck to the specimen and when the force is close to 0 N, an extra force has to be generated to break the bonding between the pins and the specimen. Due to the fact that the lower force doesn’t equal zero ΔKI is not exactly equal to KI but KIupperKIlower. But KIlower is in the range of 2 MPa*√m and therefore neglected. The upper boundary of the force is between 10 and 15 kN. It’s no problem to change the maximum force during an experiment because metal isn’t history dependent. Higher forces are also applied but the crack propagates fast which might be caused by a too large plastic zone. If the force is high enough, the plastic zone can be even seen with the micro camera (see Figure E, F and G in Appendix G). Though the specimen starts to dissipate energy substantially at a frequency of hundreds of Hertz and the maximum frequency of the draw bench is 10 Hertz, a frequency of 4 Hertz is applied. A higher frequency isn’t used because the oil in the hydraulic actuator of the draw bench once overheated. Due to a low force, specimens may need 150 000 cycles for they break. Also special clamps are developed for the fixation of the specimens with two pins (see Figure C and D in Appendix F). The micro camera and light source are of the brand Volpi and the other supporting hardware is made by Nikon (see Figure B, C and D in Appendix F). The camera is mounted on a special designed support and is placed in front of the sample. The last part of the support on which the camera is mounted to, can be moved. This enables the user to move it in 3 directions. In this way it can follow the crack tip and it can be focused. It’s not possible to rotate the camera. The micro camera itself can also be focused but this is fixed so that every image has the same scale. If it’s in focus, the software automatically scales the image so that desired lengths can be determined by the user. In our experimental set up the actual crack tip position can be located within ± 2.5 µm.. Compared with commercially used techniques this is not bad at all, because they mostly reach accuracies of decades of microns. At the beginning of an experiment the camera is placed at the notch and brought into focus. The upper and lower boundary of the force and the cycling frequency are set and the experiment can commence. After several thousands of cycles the experiment is stopped and there is looked if the crack propagated. If this is the case, a picture is made with the micro camera and the crack tip is labeled with the matching number of cycles. Subsequently, the camera is moved so that the crack tip and just a small part of the crack are imaged and it is brought into focus again. After that, the experiment is restarted again and the same procedure is repeated until the crack is propagated through the whole specimen. The procedure can be found in Appendix G, Figure A, B and C. The whole crack path is pictured in Figure D in Appendix G. 5.3 Problems with draw bench and micro camera A number of problems with the draw bench occurred during the experiments because it’s an old machine. The experiments are already time consuming because of a high number of cycles and the occurred problems made the experiments last even longer. The first problem that occurred was the overheating of the oil in the generator. This was probably caused by a too high applied frequency of the loading. After this, the frequency in the next experiments was set on 4 Hz and the problem was solved. Another problem was a human mistake. Once, the emergency button was pushed by accident. Not a big problem one should expect but when the experiment was restarted an error appeared in the control unit of the draw bench. The error could not be solved

17

right away and electricians had to solve it the next days. The above mentioned problems are fortunately not destructive for the specimens. The last problem that occurred with the draw bench was also not destructive but it was responsible for at least one wrong measurement. During one experiment the applied force was set on 15 kN but in one cycle the applied force raise to over 100 kN. There is no logical explanation for this, except that it’s an old machine and it has it own will once in a while. Fortunately this overload happened only once. The consequence of the overload is that at the crack tip a not negligible plastic zone occurs. This can be seen in the Figure 5.3.1. In this plastic zone a compressive stress will occur. Due to this, it takes much more cycles to propagate the crack than it normally does and at least one measurement is incorrect. Dependent of the size of the occurred plastic zone it’s possible that more measurements are incorrect.

Figure 5.3.1: In the upper left picture a plate with crack and applied stress can be seen. In the upper right picture the applied stress as a function of the time. If the applied stress is high, a non negligible plastic zone occurs at the crack tip (position A). In the rest of the plate the stress doesn’t reach the yield stress (position B). In the lower left picture the stress-strain response for position A and B can be found. The stress in position B remains elastic and if the plate is unloaded the stress will go back to zero. In position A the deformation is plastic and when the plate is unloaded the stress will become compressive next to the crack tip. Because the plate has to be in equilibrium when it’s unloaded a tensile stress will occur next to the plastic zone.

Next to the problems with the draw bench, also two problems with the micro camera occurred. The first one isn’t a big problem because it only caused some delay. It concerned some difficulties with the network connection between the hardware of the micro camera and the computer. The second problem was a problem with the scaling in the software. Because the scale was changed during the first experiment the crack length can’t be evaluated accurately and the results of the first experiment are useless.

18

6. Results

In chapter 3 was already proven on the hand of a 3D FE-model and an analytic formula that the 2D plane strain FE-model correctly calculates the J-integral, and thus the first stress intensity factor. This was only done for one case, so it’s not yet known if the trends of the crack path and force in the J-integral are predicted properly with the 2D FE-model. To investigate this, a plot is made of the J-integral as a function of the crack path and the force with the 2D FE-model and the analytic formula (see Figure 6.1 and 6.2). This was easily done with the Python script.

Figure 6.1: J-integral as a function of the force and the crack length for the 2D plane strain FE-model.

Figure 6.2 J-integral as a function of the force and the crack length for the analytic model for plane strain geometries.

19

As can be seen in the pictures above, the trends predicted with the 2D plane strain FEmodel are the same as those predicted with the analytic model. The J-integral is in both cases a quadratic function of the force. This can be explained by Formula 3.2 and 3.3. From Formula 3.3 it’s clear that the first stress intensity factor is a linear function of the force and in Formula 3.2 is stated that the J-integral is a quadratic function of the first stress intensity factor and thus of the force. From the same point of view can be explained that the J-integral is a polynomial of the crack length. In Figure 6.3 the error of the two J-integrals is plotted. The error is everywhere in the plot more or less the same (around 10%). So the trends are also correctly predicted with the 2D plane strain model. One may assume that the 2D model predicts the J-integral more accurately than the analytical model because the analytical model only takes into account the applied force, thickness, crack length and width of the specimen. This means that it ignores crack height and the actual inner and outer geometries.

Figure 6.3: The error of the J-integral predicted with the 2D FE-model and the one predicted with the analytic formula.

To get an idea what the displacement of the draw bench will be and how the change of the displacement is dependent of the crack length and the force, a plot is made of the crack opening displacement as a function of the force and the crack length (see Figure 6.4). It’s not possible to calculate this with the analytic model because there is no formula for derived. So a comparison between the FE-model and the analytic model isn’t possible. It can be seen that the crack opening is a linear function of the force and a polynomial of the crack length. The crack opening in the experiments will be in the order of magnitude of decades and hundreds of microns.

20

Figure 6.4: Crack opening as a function of the force and the crack length, calculated with the 2D FEmodel.

Now we definitely know that the 2D FE-model for plane strain geometries predicts the J-integral correctly, it can be used for the evaluation of the stress intensity factor in the experiments. In Figure 6.5 and 6.6 all the measurements are pictured of the specimens with a thickness of 10mm. It may be clear that all the points in Figure 6.5 are a combination of the experiments (da/dN follows from these) and the numerical analyses made by the Python script (K follows from these). The measurements are at first sight satisfying because the points are on a line in Figure 6.5 and in Figure 6.6 an exponential relation can be seen as in the theoretical results in Appendix B. The results are especially satisfying from the point of view that a number of problems occurred during the experiments as is written in chapter 5. In the results of experiments 2 and 3 a change of exponential relation can be seen in Figure 6.6 because the applied force is changed during the experiment. A least mean square fit is made on the measurements in Figure 6.5. This fit matches Paris’ law (see Figure 6.5). A good parameter for the quality of the least mean square is the variance which equals 0.0391 m/cycle. If we evaluate the unknown parameters in Paris’ law, m equals 3.371 and C equals 8.312 * 10-12 m/(MPa*√m)3.371. These results match data from literature for steel. If we look carefully to Figure 6.5 we can see that there are three measuring points of experiment 2 quite far from the fit. These are probably errors and therefore they are canceled. The newly obtained results can be found in Figure 6.7. The variance is 0.0194 m/cycle which is much better than the previous fit. The new values for m and C are 3.229 and 0.125 * 10-11 m/(MPa**√m)3.229. These results still match the values from literature for steel. These parameters may be dependent on the geometry and frequency, but because there are no other specimens and experimental setups available it’s not possible to investigate if the obtained values for C and m are the same in different experiments. To make sure that the results are really located in domain II in which Paris law is valid, the threshold and critical value for the stress intensity factor are also plotted in Figure 6.7. The threshold value can be determined with experiments but for this, an

21

extremely high number of cycles is necessary and because the maximum allowable frequency of the draw bench is only 4 Hz, this will take too long. Values for Kth for steel from literature are 6.6 MPa*√m. The critical value for the stress intensity factor can be determined with our equipment. This procedure is done after a normal experiment when the specimen is still in one piece, only the crack has propagated. When the normal experiment is stopped, the force is slowly increased (no cycling load) until the crack starts to grow. The force at this moment (52 kN) should be used in the numerical analysis to determine the critical stress intensity factor (151.62 MPa*√m). As one can see in Figure 6.7, the results are almost exactly between the Kth and Kc. This means that they are located in the regime where Paris’ law is valid.

Figure 6.5: Log-log plot of the crack propagation per number of cycles as a function of the stress intensity factor with all the measurements of the samples with a thickness of 10mm and a least mean square fit which matches Paris’ law.

22

Figure 6.6: Plot of the crack length as a function of the number of cycles with all the measurements of the samples with a thickness of 10mm.

Figure 6.7: Log-log plot of the crack propagation per number of cycles as a function of the stress intensity factor of the measurements of the first experiment with a least mean square fit that matches Paris’ law. The pink and red lines determine respectively the threshold and the critical value of the stress intensity factor.

23

Now Paris’ law is known, it is been implemented in the user subroutine. Two experiments with specimen B are done to check the user subroutine (see Figure B in Appendix A). In this way it can be checked if it predicts the correct crack path as well as the correct number of cycles. In Figure 6.8 the results of the crack paths can be found. The results of experiment 1 have to be ignored because during this experiment the emergency button was pushed by accident. It’s not known what happened with the loading at that moment but the results of the crack propagation per number of cycles suggest that this experiment has to be canceled. Therefore, only experiment 2 has to be taken into account. The simulation of experiment 2 matches the crack path very properly. A deviation can be found only at the end. This is probably due to a large plastic zone in the experiment which of course is not taken into account in the simulation because it only uses linear elastic material behavior. The results of the number of cycles as a function of the crack propagation of experiment 2 are pictured in Figure 6.9. One can clearly see the moment the load was changed from 7.5 to 6 kN. The results of the numerical simulation can be found in blue. The red line matches the results of the numerical simulations including T-stress correction which is not taken into account in this report. Therefore the red line can be ignored. Unfortunately this numerical simulation is done with different values for m and C than are determined in this report because at the time the simulations were done, we didn’t have so many data to establish Paris’ law. At that time the values of m and C were determined to be 2.98 and 0.250 * 10-11 m/(MPa**√m)2.98. The blue line matches the green line reasonable, but for the longer prediction it clearly differs more and more. This is probably because the values of m and C were not correct at that time because they were determined from only one experiment. With the new values the numerical simulation would probably match the experimental results better. The newly obtained value of C is half the value for which the simulation was done with. This means that the new numerical results will be even higher than the old ones and the slope will also be twice as high for low values of K. On the other hand, the newly obtained value of m is higher than the old one which means that slope of the new line will eventually become lower than this old one. This is exactly like in the experiment. Unfortunately it’s not possible to evaluate if the user subroutine works correctly and if the correct values of Paris’ law are established. It seems that the user subroutine on its own works fine, though there is only one prediction made with it. There are simply not enough experiments done to determine if the values in Paris law are correct.

24

Figure 6.8: Crack path to the circular hole. The results of experiment 1 can be considered as incorrect due to problems with the draw bench during experiment 1. The results of experiment 2 match the numerical simulation of experiment 2. During experiment 1 the load was set on 8 kN and at a crack length of 2200 µm. changed to 10 kN. The load in experiment 2 was set on 7.5 kN and changed to 6 kN at a crack length of 3000 µm.

Figure 6.9: Number of cycles as a function of the crack length for experiment 2 in Figure 6.8. The change of load from 7.5 kN to 6 kN can be clearly seen at a crack length of 3000 µm. The blue line is the numerical simulation without T-stress, the red line is the numerical simulation with T-stress and can be ignored in this report. The simulation matches the experiment reasonable until 5000 cycles, after that the prediction differs more and more from the experiment.

25

7. Conclusion

An experimental, numerical method is presented in this report to establish Paris’ law for steel CSN EN 100083-1+A1. Fatigue experiments with compact tension specimens in plane strain conditions are done to determine the crack propagation per number of cycles with a micro camera. After a certain number of cycles the experiments are stopped and an image is made with the micro camera. The lengths of crack propagation can easily be established with the automatically scaling in the software of the micro camera. This method allows actual crack tip positioning within ± 2.5 μm which can be rated as good compared with a number of commercially used techniques. Because it’s not possible to evaluate the first stress intensity factor accurately from analytical models of the experiments, numerical analyses were necessary. Because the geometry of a sample changes every time the crack propagates, a lot of different geometries have to be modeled and this will consume a lot of time. Therefore, a Python script was written in Complete Abacus Environment that automatically determines the J-integral for any given 2D geometry with plane strain conditions. The first stress intensity factor can easily be calculated from the J-integral for loading in mode I. The 2D numerical analyses are checked with 3D numerical analyses and analytically calculated J-integrals and it was proven that they predict J-integrals correctly. Subsequently, the experimentally determined crack propagation per number of cycles are combined with the numerically calculated first stress intensity factors to determine the parameters m and C in Paris’ law. The results are satisfying because the measuring points are almost on a straight line without any divergent values. The values of m and C are respectively 3.229 and 0.125 * 10-11 m/(MPa**√m)3.229, which match data from literature for steel. To establish if the measuring points are within the power law regime where Paris’ law is valid, the critical stress intensity factor is determined by an experiment and the threshold stress intensity factor is taken from literature. The results point out that the measurements are far enough from the critical and threshold value which means that all the measuring points are in the power law regime. Though the results seem to be satisfying, not enough experiments are done to say that the values for C and m are correct. Next to that, the parameters C and m are often load, frequency and geometry dependent so it’s recommended to perform also different experiments than only with compact tension specimens. A good recommendation is the three point bending specimen. To check if the user subroutine works properly, two experiments are done with asymmetric geometry. In this way, not only the predicted number of cycles can be checked but also the predicted crack path. Unfortunately a mistake was made during one of the two experiments so it’s only possible to verify the user subroutine on the hand of one experiment which is far too little. If nevertheless an evaluation is made, one can conclude that the predicted crack path matches the experimentally determined crack path perfectly. A small deviation can only be found at the end of the crack path. This is probably caused by a larger plastic zone at the end in the experiment. This can’t be taken into account in the user subroutine because it uses linear elastic material behavior. It’s difficult to assess the predicted number of cycles because the numerical simulation was done other values for m and C than the final ones. If, on the other hand, the results are extrapolated with the new values it seems that the trends of the number per cycles per crack propagation will be predicted correctly. It’s recommended to do more tests for the validation of the user subroutine.

26

8. Symbols Symbol a N KI m C τ r θ KII J Γ W n t u x E ν P B W ξ L

Name Crack length Number of cycles First stress intensity factor Exponent index Coefficient index Shear stress Radius Angle Second stress intensity factor J-integral Path Energy Normal vector Stress component Deformation Position coordinate in first principal direction Young’s modulus Poisson’s ratio Force Thickness Width Coordinate in first principal direction Length

Unit m MPa*√m m/(MPa*√m)m MPa m ° MPa*√m N/m M J MPa m m MPa N m m m m

27

9. Literature

[1]

Broek, D, Elementary engineering fracture mechanics, Martinus Nijhoff Publishers, Dordrecht, 1982. p. 67-90, 170-184

[2]

Nair, S.V., Tien, J.K., Bates, R.C., Buck. O., Fracture Mechanics: microstructure and micromechanisms, ASM Int., Ohio, 1989. p.111-130, 255280.

[3]

Schreurs, P.J.G., Educational material for course fracture mechanics, Eindhoven.

[4]

Roylance, D., Fatigue, Massachusetts Institute of Technology, Cambridge, 2001.

[5]

Kuzelka, J., Numerický model šíření trhliny v podmínkách cyklického zatěžování s využitím dvouparametrové lomové mechaniky, Prague, 2007.

[6]

Manuals for Abacus v6.5-5.

28

10. Acknowledgements

At first, I would like to thank Miroslav Španiel that he made my internship at the CTU in Prague possible. I also want to thank him for his supervision during the research. Secondly, I want to thank Josef Jurenka for helping me with problems that occurred during this project and for his guidance during experiments.

29

11. Appendices

30

11.1 Appendix A: Geometries of specimens

Figure A: Specimen A

31

Figure B: Specimen B

32

11.2 Appendix B: Theoretical results for fatigue crack propagation

Figure A: Log-log plot of the crack growth per number of cycles as function of the stress intensity factor. Domain I is characterized by the initiation of crack, domain II by slow, stable crack growth and domain III by instable crack growth towards failure. Paris’ law is valid in domain II and its slope is m.

Figure B: Plot of the crack length as function of the number of cyles. Domain I is characterized by the initiation of crack, domain II by slow, stable crack growth and domain III by instable crack growth towards failure.

33

11.3 Appendix C: 3D & 2D FE-model of specimen in Abacus

Figure A: 3D model

Figure B: 2D model

Figure C: Detail of mesh: rosette crack tip

34

11.4 Appendix D: Python script for fast computation of J-integrals for different geometries #Material properties: E-modulus in MegaPascal E_modulus=209e3 Poisons_ratio=0.3 #give coordinates of all outer points (clockwise) in micrometers points=[(-6.5e3,-1),(-3e3,0),(-6.5e3,1),(-9.0981e3,1.5e3),(31.5e3,1.5e3),(-31.5e3,24652.59),(-28284.3,28284.3),(15307.3,36955.2),(0.0,40e3),(15307.3,36955.2),(28284.3,28284.3),(3695 5.2,15307.3),(40e3,0.0),(36955.2,-15307.3),(28284.3,28282.43),(15307.3,-36955.2),(0.0,-40e3),(-15307.3,-36955.2),(28284.3,-28284.3),(-31.5e3,-24652.59),(-31.5e3,-1.5e3),(-9.0981e3,1.5e3)] #points=[(-10e3,-40e3),(-10e3,39.99e3),(-1.6e3,39.995e3),(10e3,40e3),(-30e3,40e3),(30e3,50e3),(30e3,50e3),(30e3,40e3),(10e3,40e3),(1.5e3,39.995e3),(10e3 ,39.99e3),(10e3,-40e3),(30e3,-40e3),(30e3,-50e3),(-30e3,-50e3),(30e3,-40e3)] #points=[ (-20e3, 40e3), (20e3, 40e3), (20e3, 10e3), (10e3, 0.0), (20e3, -20e3), (-20e3, -20e3), (-10e3, 0.0), (-10e3, 0.1e3) ] #Determine the crack tips cracktipnumbers=[1] #Determine for every cracktip the vector normal to the crackplane by defining the startpoint and endpoint of the vector: CRACKVECTOR=[[(startpointcoordx1,startpointcoordy1,startpointcoordz1) ,(endpointcoordx1,endpointcoordy1,endpointcoordz1)],[(startpointcoord x2,startpointcoordy2,startpointcoordz2),(endpointcoordx2,endpointcoor dy2,endpointcoordz2)]] CRACKVECTOR=[[(0,0,0),(0,1,0)]] #Determine for every cracktip the number of contours Number_of_contours=[5] #Determine edges to prescribe boundary conditions and forces to by giving the 2 points between an edge bc_edge_numbers=[(1,2),(8,9)] #bc_edge_numbers=[(5,6),(13,14)] bc_edge_numbers=[] #Boundary Conditions for edges: #displacement and rotation fixed: 'SET', if not fixed: 'UNSET' #Ux=displacement in x direction, Uy=displacement in y direction, URz= rotation around z axis #Fill everything in like [[Ux1, Uy1, URz1],[Ux2, Uy2, URz2]] #BC_e=[['SET','UNSET','SET'],['SET','SET','SET']] BC_e=[] #Forces on edges: #4 options: 0) no forces prescribed for this circle # 1) x force prescribed for this circle # 2) y force prescribed for this circle # 3) x and y force prescribed for this circle #Fill always [[Option1, Fx1, Fy1], [Option2, Fx2, Fy2]] in microNewton #If there is no force prescribed just fill in number 0 #F_e=[[2,0,1000e6],[0,0,0]]

35

F_e=[] #thickness of geometry in micrometers thick=10e3 #if innergeometry contains circles, give coordinates of center and 1 point on circle inner_points_c=[ [center1, point],[center2, point] ] #inner_points_c=[[(-5e3,5e3),(-4e3,5e3)],[(5e3,-5e3),(0.5e3,-4e3)]] inner_points_c=[[(-19e3,13.75e3),(-25.5e3,13.75e3)],[(-19e3,13.75e3),(-25.5e3,-13.75e3)]] #inner_points_c=[] #Boundary Conditions for circles: #displacement and rotation fixed: 'SET', if not fixed: 'UNSET' #Ux=displacement in x direction, Uy=displacement in y direction, URz= rotation around z axis #Fill everything in like [[Ux1, Uy1, URz1],[Ux2, Uy2, URz2]] BC_c=[['SET', 'UNSET', 'UNSET'],['SET', 'SET', 'UNSET']] #BC_c=[] #Forces on circles: #4 options: 0) no forces prescribed for this circle # 1) x force prescribed for this circle # 2) y force prescribed for this circle # 3) x and y force prescribed for this circle #Fill always [[Option1, Fx1, Fy1], [Option2, Fx2, Fy2] #If there is no force prescribed just fill in number 0 F_c=[[2,0,11000e6],[0,0,0]] #F_c=[] #if innergeometry contains rectangles, give coordinates of corners clockwise inner_points_r=[[corner1_1, corner2_1, corner3_1, corner4_1],[corner1_2, corner2_2, corner3_2, corner4_2]] #inner_points_r=[[(-5e3, 3e3),(5e3,3e3),(5e3,1e3),(-5e3,1e3)]] inner_points_r=[] #Boundary Conditions for rectangles: #displacement and rotation fixed: SET, if not fixed: UNSET #Ux=displacement in x direction, Uy=displacement in y direction, URz= rotation around z axis #Fill everything in like [[Ux1, Uy1, URz1],[Ux2, Uy2, URz2]] #BC_r=[['SET','SET','SET']] BC_r=[] #Forces on rectangles: #4 options: 0) no forces prescribed for this rectangle # 1) x force prescribed for this rectangle # 2) y force prescribed for this rectangle # 3) x and y force prescribed for this rectangle #Fill always [[Option1, Fx1, Fy1], [Option2, Fx2, Fy2]] #If there is no force prescribed just fill in number 0 #F_r=[[0, 0, 0]] F_r=[] #determine mesh size meshsize=3.5e3 #--------------------------------------------------------------------from part import *

36

from from from from from from from from from from

material import * section import * assembly import * step import * interaction import * load import * mesh import * job import * sketch import * visualization import *

#--------------------------------------------------------------------myModel = mdb.Model(name='SpecimenA') x=[] #--------------------------------------------------------------------mySketch = myModel.Sketch(name='Sketch Specimen', sheetSize=0.1) #Outer geometry startangles=[] endangles=[] Radii=[] new_pointsA=[] new_pointsB=[] P_radii=[] partitionpointsA=[] partitionpointsB=[] #creating straight lines of outer geometry lines=[] for i in range(len(points)-1): lines.append(mdb.models['SpecimenA'].sketches['Sketch Specimen'].Line(point1=points[i], point2=points[i+1])) lines.append(mdb.model['SpecimenA'].sketches['Sketch Specimen'].Line(point1=points[len(points)-1], point2=points[0])) #deleting lines next to cracktips selectedDeletedLines=[] for i in range(len(cracktipnumbers)): selectedDeletedLines.append(lines[cracktipnumbers[i]-1]) selectedDeletedLines.append(lines[cracktipnumbers[i]]) mdb.model['SpecimenA'].sketches['Sketch Specimen'].delete(selectedDeletedLines) #creation of arcs near cracktips and partioncircles near cracktips k=0 for i in range(len(cracktipnumbers)-1): LEN2=len(cracktipnumbers)-1-i for j in range(LEN2): DIS=((points[cracktipnumbers[i]][0]points[cracktipnumbers[j+1+i]][0])**2+(points[cracktipnumbers[i]][1]points[cracktipnumbers[j+1+i]][1])**2)**0.5 if DIS= 0: if (points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1]) >= 0: startangle=asin(abs(points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1])/radius1) else: startangle=2*pi-asin(abs(points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1])/radius1) else: if (points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1]) >= 0: startangle=pi-asin(abs(points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1])/radius1) else: startangle=pi+asin(abs(points[cracktipnumbers[i]-1][1]points[cracktipnumbers[i]][1])/radius1) if (end[0]-points[cracktipnumbers[i]][0]) >= 0: if (end[1]-points[cracktipnumbers[i]][1]) >= 0: endangle=asin(abs(end[1]points[cracktipnumbers[i]][1])/radius2) else: endangle=2*pi-asin(abs(end[1]points[cracktipnumbers[i]][1])/radius2) else: if (end[1]-points[cracktipnumbers[i]][1]) >= 0: endangle=pi-asin(abs(end[1]points[cracktipnumbers[i]][1])/radius2) else: endangle=pi+asin(abs(end[1]points[cracktipnumbers[i]][1])/radius2) startangles.append(startangle) endangles.append(endangle)

38

#defining radius (from cracktipnumber i) of the new points that need to be created for arc if ( (abs(cracktipnumbers[i]-cracktipnumbers[i-1])!=1 and abs(cracktipnumbers[i]-cracktipnumbers[endindex])!=1) and (cracktipnumbers[i]!=0 and cracktipnumbers[i-1]!=(len(points)-1))) or (cracktipnumbers[i]==0 and cracktipnumbers[i-1]!=(len(points)-1)) or (cracktipnumbers[i]==0 and len(cracktipnumbers)==1): if radius1>=1.5e3 and radius2>=1.5e3: Radius=0.01e3 P_radius=1.5e3 else: Radius=min([radius1,radius2])/150 P_radius=min([radius1,radius2]) if (abs(cracktipnumbers[i]-cracktipnumbers[i-1])==1 and abs(cracktipnumbers[i]-cracktipnumbers[endindex])!=1) or (cracktipnumbers[i]==0 and cracktipnumbers[i-1]==(len(points)-1) and abs(cracktipnumbers[i]-cracktipnumbers[endindex])!=1): if radius1>=3e3 and radius2>=1.5e3: Radius=0.01e3 P_radius=1.5e3 else: Radius=min([radius1/2,radius2])/150 P_radius=min([radius1/2,radius2]) if ((abs(cracktipnumbers[i]-cracktipnumbers[i-1])!=1 and abs(cracktipnumbers[i]-cracktipnumbers[endindex])==1) and (cracktipnumbers[i]!=0 and cracktipnumbers[i-1]!=(len(points)-1))) or (cracktipnumbers[i]==(len(points)-1) and cracktipnumbers[endindex]==0): if radius1>=1.5e3 and radius2>=3e3: Radius=0.01e3 P_radius=1.5e3 else: Radius=min([radius1,radius2/2])/150 P_radius=min([radius1,radius2/2]) if (abs(cracktipnumbers[i]-cracktipnumbers[i-1])==1 and abs(cracktipnumbers[i]-cracktipnumbers[endindex])==1) or (cracktipnumbers[i]==0 and cracktipnumbers[i-1]==(len(points)-1) and abs(cracktipnumbers[i]-cracktipnumbers[endindex])==1): if radius1>=3e3 and radius2>=3e3: Radius=0.01e3 P_radius=1.5e3 else: Radius=min([radius1/2,radius2/2])/150 P_radius=min([radius1/2,radius2/2])

Radii.append(Radius) P_radii.append(P_radius) #defining points for creation of arc new_pointA=(points[cracktipnumbers[i]][0]+Radius*cos(startangle), points[cracktipnumbers[i]][1]+Radius*sin(startangle)) new_pointB=(points[cracktipnumbers[i]][0]+Radius*cos(endangle), points[cracktipnumbers[i]][1]+Radius*sin(endangle)) new_pointsA.append(new_pointA) new_pointsB.append(new_pointB)

39

partitionpointA=(points[cracktipnumbers[i]][0]+P_radius*cos(startangl e), points[cracktipnumbers[i]][1]+P_radius*sin(startangle)) partitionpointB=(points[cracktipnumbers[i]][0]+P_radius*cos(endangle) , points[cracktipnumbers[i]][1]+P_radius*sin(endangle)) partitionpointsA.append(partitionpointA) partitionpointsB.append(partitionpointB) #creation of arc mdb.model['SpecimenA'].sketches['Sketch Specimen'].ArcByCenterEnds(points[cracktipnumbers[i]], new_pointB, new_pointA)

if k==1: getWarningReply('2 or more cracktips are to near to each other, try other geometry', (YES,NO)) mdb.close()

else: for i in range(len(cracktipnumbers)): if cracktipnumbers[i]==(len(points)-1): end=points[0] else: end=points[(cracktipnumbers[i]+1)] if (i+1) == len(cracktipnumbers): endindex=0 else: endindex=i+1 if (abs(cracktipnumbers[i]-cracktipnumbers[i-1])==1 and cracktipnumbers[i-1]!=cracktipnumbers[len(cracktipnumbers)-1]) or (cracktipnumbers[i]==0 and cracktipnumbers[i-1]==(len(points)-1)): mdb.model['SpecimenA'].sketches['Sketch Specimen'].Line(new_pointsA[i], new_pointsB[i-1]) else: mdb.model['SpecimenA'].sketches['Sketch Specimen'].Line(new_pointsA[i], points[(cracktipnumbers[i]-1)]) if (abs(cracktipnumbers[i]-cracktipnumbers[endindex])==1 and cracktipnumbers[endindex]!=cracktipnumbers[0]) or (cracktipnumbers[endindex]==0 and cracktipnumbers[i]==(len(points)1)): x.append(1) else: mdb.model['SpecimenA'].sketches['Sketch Specimen'].Line(new_pointsB[i], end)

#Inner geometry for i in range(len(inner_points_c)): mdb.model['SpecimenA'].sketches['Sketch Specimen'].CircleByCenterPerimeter(inner_points_c[i][0], inner_points_c[i][1])

40

for i in range(len(inner_points_r)): for j in range(len(inner_points_r[i])): mdb.model['SpecimenA'].sketches['Sketch Specimen'].Line(inner_points_r[i][j-1],inner_points_r[i][j])

#--------------------------------------------------------------------#Creation of part myPart = myModel.Part(name='Part A', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY) mdb.models['SpecimenA'].parts['Part A'].BaseShell(sketch=mySketch)

#--------------------------------------------------------------------#Creation of partition partitionSketch=mdb.models['SpecimenA'].Sketch(name='partition1', sheetSize=0.1e6) for i in range(len(cracktipnumbers)): mdb.models['SpecimenA'].sketches['partition1'].ArcByCenterEnds(points [cracktipnumbers[i]], partitionpointsB[i], partitionpointsA[i])

CrackRegionFace=mdb.models['SpecimenA'].parts['Part A'].faces.findAt((partitionpointsA[0][0], partitionpointsA[0][1], 0.000000),) mdb.models['SpecimenA'].parts['Part A'].PartitionFaceBySketch(faces=CrackRegionFace,sketch=partitionSketc h)

#--------------------------------------------------------------------# assembly import assembly mdb.models['SpecimenA'].rootAssembly.Instance(name='Assembly A', part=mdb.models['SpecimenA'].parts['Part A']) #--------------------------------------------------------------------#Creation of mesh import mesh

model_edges=mdb.models['SpecimenA'].rootAssembly.instances['Assembly A'].edges model_points=mdb.models['SpecimenA'].rootAssembly.instances['Assembly A'].vertices

41

InstancesA=mdb.models['SpecimenA'].rootAssembly.instances['Assembly A']

mdb.models['SpecimenA'].rootAssembly.seedPartInstance(regions=(Instan cesA,), size=meshsize) elemType1 = mesh.ElemType(elemCode=CPE8, elemLibrary=STANDARD) elemType2 = mesh.ElemType(elemCode=CPE6M, elemLibrary=STANDARD) crack_regions=[] for i in range(len(cracktipnumbers)): ed11=model_edges.findAt( (((partitionpointsA[i][0]new_pointsA[i][0])/2+new_pointsA[i][0]), ((partitionpointsA[i][1]new_pointsA[i][1])/2+new_pointsA[i][1]), 0.0000000), ).index ed21=model_edges.findAt( (((partitionpointsB[i][0]new_pointsB[i][0])/2+new_pointsB[i][0]), ((partitionpointsB[i][1]new_pointsB[i][1])/2+new_pointsB[i][1]), 0.0000000), ).index if endangles[i]