Implementation of Special Singular Finite Elements for Cracks along Adhesive Interfaces in Mode III

28 Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Ingeniería Aeroespacial Implementation of Special Singular Finite Elements for Cracks along Adhesive Interfaces in Mode III Autor: Antonio Manuel Alonso Arias Tutores: Vladislav Mantič Leščišin, Luis Arístides Távara Mendoza Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Grado Ingeniería Aeroespacial Implementation of Special Singular Finite Elements for Cracks along Adhesive Interfaces in Mode III Autor: Antonio Manuel Alonso Arias Tutores: Vladislav Mantič Leščišin, Luis Arístides Távara Mendoza Profesores Titulares Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Carrera: Implementation of Special Singular Finite Elements for Cracks along Adhesive Interfaces in Mode III Autor: Antonio Manuel Alonso Arias Tutores: Vladislav Mantič Leščišin, Luis Arístides Távara Mendoza El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la n de: El Secretario del Tribunal Fecha: Acknowledgements Scientific knowledge and engineering only make sense when they are bound to make people more free and happier. When technology improves human lives, it allows them to look beyond their own survival and to aim for discovering the world around them. Hence, they become free to make decisions by themselves, and ultimately prevent others from making decisions on behalf of them (and usually, against them). If technology has such a democratic purpose, why is it that elitism makes its way to many a technician’s heart? I want to thank those who encouraged me to maintain my thirst for knowledge amidst this impassible world. They understood my failures and corrected my errors because they saw potential in me and my generation. Among them are my tutors Mar, Luis, and Vladislav; and professors Miguel Pérez-Saborid, Cristina Arévalo, and Eva Pérez. A special thanks to Alberto, whose simulations were crucial to the research and to Alejandro Estefani for organising an Abaqus seminar. As Ludwig van Beethoven once said: ‘I know no other sign of superiority than Goodness’. Finally, I want to thank those who were also in darkness but made this world worth fighting for. The light shines in the darkness, and the darkness has not understood it. Antonio Manuel Alonso Arias Sevilla, 2022 I Resumen El MEF puede aplicarse a la Mecánica de la Fractura para simular tensiones en el entorno del borde de grieta; así se puede calcular el Factor de Intensificación de Tensiones. Los elementos lineales aportan tensiones constantes en su interior, por lo cual no son apropiados para captar la tendencia al infinito, requiriéndose mallas densas. Para reducir el coste computacional, se diseñan Elementos Singulares, cuyas funciones de forma y/o disposición de los nodos permite captar mejor las tensiones en torno al borde de grieta con menor número de elementos. Aun así, un comportamiento correcto de cualquier elemento singular es muy dependiente del mallado de la pieza: forma, disposición y tamaño de los elementos alrededor de la región mallada con elementos singulares, y el tamaño relativo entre elementos y grieta entre otros. El objetivo de este trabajo es desarrollar un breve manual de uso de un nuevo Elemento Singular con singularidad logarítmica adaptado para interfases adhesivas en modo III. Como objetivo secundario, el proceso de pruebas se generaliza para extrapolarse fácilmente. Finalmente, los códigos y macros generados se adjuntan para que el proceso se aplique fácilmente a nuevos Elementos Singulares. III XContents Study 1 148 Study 2 152 Study 3 161 Figures Index 165 Tables Index 169 Scripts Index 171 References 173 Notation We will be using Einstein’s notation BVP Boundary Value Problem SE Singular Element ul Unit of length up Unit of pressure KI,II,III Stress Intensity Factor for modes I, II or III KIC,IIC,IIIC Fracture toughness of the material for modes I, II or III KGeneralized Intensity Factor σi j Cauchy stress tensor Xi j Volumetric forces affecting a solid DSolid domain ∂DSolid’s contour region εi j Strain tensor EYoung Modulus σEYield strength νPoisson’s ratio KBulk modulus GShear modulus λLamé’s first parameter µLamé’s second parameter ∂DtSolid’s contour region where tension is known ∂DuSolid’s contour region where displacements are known ∂Dk Solid’s contour region where there is a relation between displacements and stresses TTraction vector LQQuarter-point element size LNNormal element size LNi size of the ith normal element layer which encompasses the other (i-1)th singular elements LTTransition element size LS Size of the singular elements (either normal or transition elements) LSRZ Size of the SRZ u,ΨDisplacement fields ξ,ηNatural coordinates of an finite element αiParameters of a crown XI 1 The Case Study If one must fight or create, it is necessary that this be preceded by the broadest possible knowledge. Karel Čapek The purpose of this work is finding out the positive traits of a mesh which utilizes the Singular Elements developed by Manuel Romero in [ 17 ] for 2D cracks along adhesive boundary layers in mode III. We will focus on the traits of the elements surrounding the crack: both the singular elements, and the non-singular elements around them (which generally form one or several crowns, e.g. [ 34 ], [ 6 ]). Finally, we will figure out guidelines of size, shape and placing (among others) of the elements, which will allow us to compute a more realistic solution with the lesser number of elements possible. Following those guidelines will reduce the operational cost of solving the problem thus reducing solving time. 1.1 Hypotheses Crack y xz 1 ul 2 ul Crack tip Figure 1.1 Plate and Reference System. Our case study has the following shape and characteristics: • Brittle linear elastic material with no plastification. This models a thin (orthotropic with principal directions tangential and normal to the surface) adhesive layer that joins two plates of other (isotropic and homogeneous) materials, for instance composite material laminates (of isotropic behavior). This model (and its FEM approach) has been proven valid in [22], [33] and [21]. • Body approximated by a 2D Shell of base 1ul and height 2ul. • Crack growing in the longitudinal half of the plate and expanding from one lateral edge to the center of the plate. • The coordinate system will be as described in the diagram on the left. • Neumann boundary conditions in the top edge (with value 1up ), bottom edge (with value −1up ), lateral edges and open crack lips (both with null value). They model a continuous load on the top and bottom edges of the plate perpendicular to its plane. • Robin boundary conditions in the crack growth direction, with the elastic constant ( k ) as a parameter. They model an adhesive interface which maintains the two parts joined at the right half but has been broken at the left half. This is equivalent to a homogeneous 1 2Chapter 1. The Case Study and constant distribution of springs. Ultimately this relates the perpendicular stress at each pair of correlative points in the upper and lower edges with their relative displacement (see Hooke’s low later). From all of the above, we can infer: • From the above we can infer a load skew-symmetry, hence we can study only the top half with minor changes. For instance, the half-model springs will have twice the stiffness than the whole-model springs. Stress on a node in Robin edges is equal to the elastic constant of the spring times its elongation (which we assume equal to the node displacement). The half-model elongation is half of the whole-model elongation since the origin of coordinates is in the undeformed plate, halfway between the separated sides of the plate. Hence, the half-model spring stiffness should be doubled to compensate this fact (if we want to calculate the real displacement field). • Finally, we can also infer that the crack grows under mode III (anti-plane problem), separating the upper and lower crack lips and sides of the Robin edges. These edges are still joined together by the adhesive layer which remains intact up until the crack growth. • Let us accept that we can define a Generalized Stress Intensity Factor K (which would not exactly be considered KIII due to the spring distribution) in order to compare it to the Stress Intensity Factor. As remarked by [ 27 ], these models have a practical use on modern aviation (e.g. F18) as most composite material structures have plies that need adhesives to be joined, and even some non-composites are also bonded. Aircrafts like B787 have carbon laminate wings attached to fiber glass leading edges, and A350 has wings mainly composed of carbon composites bonded to monolithic materials that compose the flaps. Crack Crack tip +1 Neumann boundary condition -1 Neumann boundary conditions null Neumann boundary condition Robin boundary condition Figure 1.2 Whole model boundary conditions. Crack Crack tip Neumann boundary condition null Neumann boundary condition Robin boundary condition Figure 1.3 Half model boundary conditions. As in other mode III cases, the stress tendencies in the vicinity of the crack are only dependant on the geometry of the specimen and the loads. However, unlike habitual mode III cases, in which stresses tend to infinity near the crack tip, the presence of springs make displacement (and according to Hooke’s law, stresses) finite on the crack tip. More precisely, according to [ 4 ], displacements are proportional to a series of terms of the shape rnlog(r) and rnlog2(r) , which (although indeterminate for r=0 ) has a finite limit. Its gradient, 1.2 Process 3 however, will grow proportionally to log(r) + 1 , which has a singularity for r=0 . This type of tendency can be difficult to capture in FEM with normal means (that is, no special elements) due to the logarithmic singularity in stresses. A more thorough study of those cases can be found in [15] or [18]. 1.2 Process Small changes in the code can make the following process applicable to cases with different boundary conditions. Also, for this study we are testing a Singular Element developed in [ 17 ] using the base of functions found in [ 4 ] [ 1 ]. Such element will allow us to properly fit our calculations to the tendency of displacements around the crack tip (employing a similar principle to the Quarter Point Element). The goal is to work out a Generalized Stress Intensity Factor for Mode III with Neumann-Robin conditions ( KIIIN−R ) in order to determine whether the crack advances by comparing its value to the adhesive layer’s fracture toughness. However, there are still no fixed guidelines for the element usage conditions. Firstly, we ought to determine the meshing properties that would optimize the element’s accuracy; that is to say, a handbook. More precisely, we will be considering how variations in the following parameters affect the FEM Outcome of the model (they will be further explained in the following chapter): • Shape of the elements around the crack tip: depending on the direction of the elements edges we will have worse or better results since if affects the mesh quality and the number of nodes in the approaching direction. • Number of elements around the crack tip: we will consider cases with 4, 6, 8, 12 and 24 in case a higher number of elements improves the approximation. •Size of the elements in proportion to the crack. • Effect of the use of a crown of non-singular elements surrounding the singular elements around the crack. • Effect of the use of two crowns of non-singular elements surrounding the singular elements around the crack. • Proportion between the size of singular elements and the sizes of non-singular elements composing each crown. • Proportion between the number of singular elements and the number of non-singular elements composing each crown. • Size of the singular-element region in proportion to the half-hexagon around it employing the optimal proportion among crown elements. •Size of the singular-element region in proportion to the total crack length. • Future use of singular elements (e.g., Quarter point or others) instead non-singular elements surrounding the singular elements. The result will be considered favorable if using singular elements arranged in a certain configuration improves the result obtained without using no singular elements. That is, if a coarse or fine meshing with singular elements converges faster than a finer meshing without singular elements. Since there is no analytical solution for this precise case values, we will be taking as an exact solution the result of a convergence study with fine meshing, using N=1024 isosceles elements on each side of the plate of the same size in [ 22 ], and considered acceptable since it has been obtained by several meshing types. Some codes, Abaqus macros, and procedures that may be useful have been developed. Even if they are just tools used along the process, they are nonetheless at the service of the University for further studies. in in in in in in in in in in in in sin cos length if if if 2.8 Our element 17 Niu ∂Ni ∂ξ ∂Ni ∂η i (1−η)1−ξ−ξln(ξ) ln(1/2)(1−η)−1−1 ln(1/2)−ln(ξ) ln(1/2)−1−ξ−ξln(ξ) ln(1/2)1 (1−η)2ξln(ξ) ln(1/2)2(1−η)ln(ξ) ln(1/2)+1 ln(1/2)−2ξln(ξ) ln(1/2)2 (1−η)ξ−ξln(ξ) ln(1/2)(1−η)1−ln(ξ) ln(1/2)−1 ln(1/2)−ξ−ξln(ξ) ln(1/2)3 ηξ−ξln(ξ) ln(1/2)η1−ln(ξ) ln(1/2)−1 ln(1/2) ξ−ξln(ξ) ln(1/2)4 η2ξln(ξ) ln(1/2)2ηln(ξ) ln(1/2)+1 ln(1/2) 2ξln(ξ) ln(1/2)5 η1−ξ−ξln(ξ) ln(1/2)η−1−1 ln(1/2)−ln(ξ) ln(1/2)η1−ξ−ξln(ξ) ln(1/2)6 (2.67) After obtaining the functions, they can be collapsed into a triangular element in the real space by making the coordinates of two of its nodes equal, which will affect the element’s behaviour. The process is an analogous to Barsoum’s process [6]. It can be illustrated as in figure 2.7. Figure 2.10 Collapse of the 2D square element into a 5 node triangle. As a result of the collapse of the nodes, shape functions number 1 and 6 must fulfill: [N1(ξ,η)+N6(ξ,η)]ξ−ηspace = [N1(x,y)]real space (2.68) Finally, to integrate each entry of the stiffness matrix, the process is analogue to that of turning a 4-noded element into a triangular one. The jacobian is chosen to affect the nodes needed, not affecting the quadraticbehaviored ones; and choosing to make equal the displacements of nodes 1 and 6 (as implied before). We end up with the following equation: 18 Chapter 2. Previous knowledge Ki j5x5=G      Z1 0Z1 0        ∂(N1+N6) ∂ξ ∂(N1+N6) ∂η ∂N2 ∂ξ ∂N2 ∂η . . .. . . ∂N5 ∂ξ ∂N5 ∂η       [J]−1T[J]−1 ∂(N1+N6) ∂ξ ··· ∂N5 ∂ξ ∂(N1+N6) ∂η ··· ∂N5 ∂η !|J|dξdη      5x5 (2.69) Stresses inside the element and the Robin matrix can be both obtained in a relatively similar manner. The integration was done in Wolfram’s Mathematica. And the whole result is added in the MATLAB code, in its own sub-section of the main program. Thanks to this, we can apply the Sanz-Herrera’s rendition of Alberty-Carstensen-Funken code to our analysis changing the stiffness matrix’s entries for our singular element equations. Lastly, what utterly motivated my work was the fact that the last of Romero’s studies did meet the expectations (that is, his element converged faster than CST’s) even if his elements were considerable deformed due to the meshing. This may have led to a worsening of the results. However, even with such drawback, results were improved. The question is now clear: can we make convergence faster with different kinds of meshing? That is to say, it is expected that different configurations can improve current conclusions, but firstly we need to know what are the beneficial traits of a mesh. Once they are identified, the study can be remade and presumably, it will show higher convergence rate. Figure 2.11 Meshes types employed up until now.. For 6 and 8 singular elements, rosette and transitioning elements are not specially deformed. However, as the number of singular elements grows, deformation becomes a problem. 2.9 Conditions of use As we can find in several sources, (for instance [ 16 ], [ 11 ], [ 20 ]), a beneficial mesh trait is the use of equilateralshaped elements, in our case equilateral triangles (or squares in other cases). Most elements are designed to bear a certain amount of distortion in order to fit more complex geometries, however an excessive deformation would dramatically affect the convergence of the solution and distort the final results of a mesh. Concepts like aspect ratio are used to calculate a mesh quality by certain softwares (among which we find Abaqus and ANSYS). Degenerate Elements augment the inaccuracy of FEA, and since our element has been obtained by 2.9 Conditions of use 19 collapsing an equilateral element (hence introducing a bit of distortion in its design), and the test meshes did also stretch too much those singular elements and the following transition elements, changing this parameter could substantially affect our result. This is what Romero suggested as changing the mesh type, that is to say, finding a configuration that improves the element behaviour close to the crack tip and does not distort the solution close enough to the singular elements as to change its behavior. Furthermore, some suggestions on element size and proportion appear in Abaqus handbooks and forums. It is recommended to keep the size of the non-singular elements around the quarter-point region as 0.7 times the singular-element size. The research in [ 19 ] remarks that quarter-point elements studies give a more accurate stress intensity factor than its quadrilateral counterpart, although specifies that natural isoparametric triangles work better than collapsed shapes. Moreover, remarks that its span angle should not exceed 60 , thus rendering our already studied case of 2 singular elements in doubt, although reinforcing the idea of using a rosette configuration is still desireable. Such elements are suggested to occupy from 5% to 25% of the crack’s length (which is later proven to be roughly the same for our element). Transition elements, those placed after the singular elements, are proposed to have a certain size in proportion to the quarter-point elements to correctly sense the singularity, being specially helpfull (according to some like [ 28 ]) when the proportion between the quarter-point is much smaller than the crack length. Our main difference with the later studies is that we will be comparing the solution ( K ) not to an analytic solution of the problem but to another numeric solution, obtained with a fine mesh which has already converge (with the number of elements). As stated in the same paper, a 1% accuracy would be more than reasonable. This explains the lack of error curves in the work, since we are comparing different shapes, not significantly augmenting the number of elements, so no shape is presumably better to compare. Figure 2.12 Taken from [19]. The same research also defined several zones: zones of dominance and zones of representation. The zones of dominance are: singular dominant zone (where stresses are dominated by the singularity, in the vicinity of the crack tip), non-singular dominant zone (where stresses are not relevantly affected by the singularity, far from the crack tip), and dual dominance zone (where neither component dominates, in-between the two others). The zones of representation are: singular stress component representation zone (the zone which we can model with singular elements), non-singular constant component stress component representation zone (the zone which can be modeled with non-singular, constant elements), and the non-singular linear stress component representation zone (which is understood to transition between both). Both zones of dominance and zones of representations should be balanced to obtain an optimal result. In our case, however, we will be centered in using constant elements instead of linear elements in the transition zone, supposing that linear behaviour will improve the solution and trying to achieve a good transition thanks to the element’s shape and arrangement instead of linear shape function. Ultimately, we will be studying the same regions but with a slightly different approach, and after the data is analysed and best-fitting data is chosen, we can estimate the real length of the zones with our final results (whose dimensions are summed in the last table) and the equations provided by Lim (variables described in notation): LSRZ =LQ+∑LS=LQ+∑LT+∑LN. 3 First Parametric Study: effect of one and two crowns surrounding the singular elements 3.1 Parameter definitions Let us begin by defining the parameters we will be modifying to determine what kind of meshing is optimal: 3.1.1 Singular element size Let us call a1 to the radius of the singular element: the distance between the crack tip and one of the vertices on the opposite side of the element (so that the whole element is encircled). For the sake of simplicity, let us call SER (Singular-Elements Region) to this region. α1 α1α1 Figure 3.1 Different values of a1. For the case of 6singular elements around the crack tip, with α1increasing towards the right: α1α1α1 Figure 3.2 Different values of a1for six singular elements. 21 22 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 3.1.2 Number of singular elements We will study the following cases: 4, 6, 8, 12 and 24 elements in a rosette pattern. Let us call n the number of elements in the neighborhood of the crack. We end up with these mesh options: n=6 n=4 n=24 n=8 n=12 Figure 3.3 Different values of nfor six singular elements. 3.1.3 Size of elements in the first crown Let us assume, as stated in the previous knowledge section, that using one crown of elements (that is, a circular crown surrounding the singular elements) will be beneficial. By now, it will be meshed with non-singular elements. The two circumferences that define the crown are defined by radii a1 and a2 . Thus a1 is both the radius of the outer circumference of the SER and also of the inner circumference of the first crown. In the image, a2increases while a1remains constant, although they are independent from one another. α1 α2 α2α2 α1α1 Figure 3.4 Different values of a2for six singular elements. 3.1.4 Type of the first crown type 1: type 2: type 4: type 5: type 3: Figure 3.5 Different values of a2for six singular elements. The different feasible shape of the elements in the first crown will be referred to as ’Type of the first crown’. We will study the following types of meshing, which are now introduced for n=6 but will be generalized in the following section (it is easily extrapolatable to any n). We can differentiate between radial crowns (types 1, 2, 3 and 4, in which the singular elements edges continue along the crown) and nonradial crowns (type 5, in which said edges do not continue). Also we can differentiate between simple crowns (types 1, 2 and 3, in which the number 3.1 Parameter definitions 23 of elements in the outer circumference is the same as in the inner one) and double crowns (types 4 and 5, in which the number of elements in the outer circumference is higher than those in the inner one). 3.1.5 Size of elements in the second crown Akin to the circumference of radius a2used for the first crown, a third circumference with radius a3is used to define the second crown. The radius increases to the right keeping a1 and a2 constant. It is worth noting that the study can be conducted with or without this crown (as we will see later). α3α3α3 α2α2α2 α1α1α1 Figure 3.6 Increasing values of a3. 3.1.6 Type of the second crown Akin to the first crown, a new one with the same types can be introduced between a2 and a3 . Let them be defined by the same codes (1, 2, 3, 4 and 5) meaning almost the same but with the following changes. In the second crown’s inner circumference there may be n elements (in case of first crown types 1, 2, 3), 2n elements (in case of the first crown type 4), or n+2 elements (in case of the first crown type 5). Hence, the ’type of the second crown’will have a slightly different shape depending on the first crown. Before further explanation in the following section, the changes are illustrated here. In the appendix these diagrams are shown for the rest of case studies (4,8,12 and 24 elements). 1st crown type: 1 2nd crown type: 1 1st crown type: 1 2nd crown type: 4 1st crown type: 1 2nd crown type: 2 1st crown type: 1 2nd crown type: 3 1st crown type: 1 2nd crown type: 5 Figure 3.7 Second crown types for a type-1 first crown. 30 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 12 s.unsetPrimaryObject() 13 p = mdb.models[ Model−1 ]. parts [ Part−1 ] 14 session .viewports[ Viewport: 1 ]. setValues ( displayedObject =p) 15 del mdb.models[ Model−1 ].sketches [ __profile__ ] Script 3.4 Step 01 in Python. • Step 02: Make the Half-Hexagon Cutting of the half-hexagon (from coordinates created assuming Origin of Coordinates on the crack tip, else we must add the crack tip coordinates plus the half-hexagon ones), the auxiliary lines from the upper vertices of the plates to the upper vertices of the half-hexagon, and the vertical symmetry line. Part module → Partition Face: Sketch → Select face, then draw lines. In yellow, the lines sketched joining auxiliary points. Figure 3.16 Result of the Step 02. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2f1, e1, d2 = p. faces , p.edges, p.datums 3t = p.MakeSketchTransform(sketchPlane=f1[0], sketchPlaneSide=SIDE1, origin=( 40.0, 0.5, 0.0) ) 5s = mdb.models[ Model−1 ].ConstrainedSketch(name= __profile__ , sheetSize =2.82, 6gridSpacing=0.07, transform=t) 7g, v, d, c = s.geometry, s. vertices , s.dimensions, s. constraints 8s. setPrimaryObject(option=SUPERIMPOSE) 9p = mdb.models[ Model−1 ]. parts [ Part−1 ] 10 p. projectReferencesOntoSketch (sketch=s, filter =COPLANAR_EDGES) 11 s.Line(point1=(0.0, -0.51), point2=(0.0+xpaux[0], -0.5+ypaux[0])) 12 s.Line(point1=(0.0+xpaux[0], -0.5+ypaux[0]), point2=(0.0+xpaux[1], -0.5+ypaux[1])) 13 s.Line(point1=(0.0+xpaux[1], -0.5+ypaux[1]), point2=(0.0+xpaux[2], -0.5+ypaux[2])) 14 s.Line(point1=(0.0+xpaux[2], -0.5+ypaux[2]), point2=(0.0+xpaux[3], -0.5+ypaux[3])) 15 s.Line(point1=(0.0+xpaux[3], -0.5+ypaux[3]), point2=(0.0+xpaux[4], -0.5+ypaux[4])) 16 s.Line(point1=(0.0+xpaux[4], -0.5+ypaux[4]), point2=(0.0, -0.51)) 17 s.Line(point1=(0.5, 0.5), point2=(0.0+xpaux[1], -0.5+ypaux[1])) 18 s.Line(point1=(-0.5, 0.5), point2=(0.0+xpaux[3], -0.5+ypaux[3])) 19 s.Line(point1=(0.0, -0.5), point2=(0.0, 0.5)) 20 p = mdb.models[ Model−1 ]. parts [ Part−1 ] 21 f = p. faces 22 pickedFaces = f .getSequenceFromMask(mask=( [#1 ] , ) , ) 23 e, d1 = p.edges, p.datums 3.3 Mesh generation 31 24 p. PartitionFaceBySketch ( faces=pickedFaces, sketch=s) 25 s.unsetPrimaryObject () 26 del mdb.models[ Model−1 ].sketches [ __profile__ ] Script 3.5 Step 02 in Python. • Step 03: Make Circumferences Cutting of the two circumferences centered in the crack tip. The inside region corresponds to singular elements, and the in-between region to the first crown. If we were to add a second crown,anew exterior circumference should be cut. Part module → Partition Face: Sketch → Select face (inside half hexagon), then draw two/three circumferences centered on the midpoint of the shell’s lower edge and with radii α1 , α2 , α3 times the hexagon side-length.. Do not use the third command when there is only one crown. Figure 3.17 Result of Step 03. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2f, e1, d2 = p. faces , p.edges, p.datums 3t = p.MakeSketchTransform(sketchPlane=f[1], sketchPlaneSide=SIDE1, origin=( 40.097222, 0.096225, 0.0) ) 5s1 = mdb.models[ Model−1 ].ConstrainedSketch(name= __profile__ , 6sheetSize =1.11, gridSpacing=0.02, transform=t) 7g, v, d, c = s1.geometry, s1. vertices , s1.dimensions, s1. constraints 8s1. setPrimaryObject(option=SUPERIMPOSE) 9p = mdb.models[ Model−1 ]. parts [ Part−1 ] 10 p. projectReferencesOntoSketch (sketch=s1, filter =COPLANAR_EDGES) 11 s1.CircleByCenterPerimeter(center=(x_ct, y_ct), point1=(x_ct+rci, y_ct)) 12 s1.CircleByCenterPerimeter(center=(x_ct, y_ct), point1=(x_ct+rco, y_ct)) 13 s1.CircleByCenterPerimeter(center=(x_ct, y_ct), point1=(x_ct+rco2, y_ct)) 14 p = mdb.models[ Model−1 ]. parts [ Part−1 ] 15 f = p. faces 16 pickedFaces = f.getSequenceFromMask(mask=(’[#a ]’, ), ) 17 e, d1 = p.edges, p.datums 18 p. PartitionFaceBySketch ( faces=pickedFaces, sketch=s1) 19 s1.unsetPrimaryObject() 20 del mdb.models[ Model−1 ].sketches [ __profile__ ] Script 3.6 Step 03 in Python. 32 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements • Step 04: Remove the SER. Removal of the inner circle, which will be substituted byadistribution of singular elements later, when the case is simulated in MATLAB. Part Module → Remove Faces, then select both halves of the innermost circle. Figure 3.18 Result of Step 04. 1session .viewports[ Viewport: 1 ]. view. setValues (nearPlane=2.63361, 2farPlane=3.02325, width=1.10348, height=0.554118, cameraPosition=( 30.0496831, 0.207573, 2.82843), cameraTarget=(0.0496831, 0.207573, 0)) 4p = mdb.models[ Model−1 ]. parts [ Part−1 ] 5f1 = p. faces 6p.RemoveFaces(faceList = f1[5:6]+f1[7:8], deleteCells=False) Script 3.7 Step 04 in Python. • Step 05: Cut the first crown. Cutting of the first crown shape according to its type. Also, cutting of the second crown shape if necessary. The figure cutting uses straight lines to replicate the edges of the elements as imposed by the crown type. This can be easily achieved since the edges of said elements coincide with segments that connect the vertices of regular polygons inscribed inside the pre-existing circumferences. Hence, we can just calculate the vertices coordinates and select them as the starting/ending points of segments (according to the desired crown type). Part Module → Select faces between remaining circumferences, then draw segments from coordinates. 3.3 Mesh generation 33 Figure 3.19 Close-up of the resulting crown. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2f, e1, d2 = p. faces , p.edges, p.datums 3t = p.MakeSketchTransform(sketchPlane=f[1], sketchPlaneSide=SIDE1, origin=( 4y_ct, 0.082041, 0.0) ) 5s = mdb.models[ Model−1 ].ConstrainedSketch(name= __profile__ , sheetSize =0.91, 6gridSpacing=0.02, transform=t) 7g, v, d, c = s.geometry, s. vertices , s.dimensions, s. constraints 8s. setPrimaryObject(option=SUPERIMPOSE) 9p = mdb.models[ Model−1 ]. parts [ Part−1 ] 10 p. projectReferencesOntoSketch (sketch=s, filter =COPLANAR_EDGES) 11 12 # First crown 13 s.Line(point1=(x_ct+x2[0] ,y_ct+y2[0] ), point2=(x_ct+ x1[1],y_ct+ y1[1])) 14 s.Line(point1=(x_ct+x1[1] ,y_ct+y1[1] ), point2=(x_ct+ x2[1],y_ct+ y2[1])) 15 s.Line(point1=(x_ct+x2[1] ,y_ct+y2[1] ), point2=(x_ct+ x1[2],y_ct+ y1[2])) 16 s.Line(point1=(x_ct+x1[2] ,y_ct+y1[2] ), point2=(x_ct+ x2[2],y_ct+ y2[2])) 17 s.Line(point1=(x_ct+x2[2] ,y_ct+y2[2] ), point2=(x_ct+ x1[3],y_ct+ y1[3])) 18 s.Line(point1=(x_ct+x1[3] ,y_ct+y1[3] ), point2=(x_ct+ x2[3],y_ct+ y2[3])) 19 s.Line(point1=(x_ct+x1[3] ,y_ct+y1[3] ), point2=(x_ct+ x2[4],y_ct+ y2[4])) 20 s.Line(point1=(x_ct+x2[4] ,y_ct+y2[4] ), point2=(x_ct+ x1[4],y_ct+ y1[4])) 21 s.Line(point1=(x_ct+x1[4] ,y_ct+y1[4] ), point2=(x_ct+ x2[5],y_ct+ y2[5])) 22 s.Line(point1=(x_ct+x2[5] ,y_ct+y2[5] ), point2=(x_ct+ x1[5],y_ct+ y1[5])) 23 s.Line(point1=(x_ct+x1[5] ,y_ct+y1[5] ), point2=(x_ct+ x2[6],y_ct+ y2[6])) 24 # Second crown 25 s1.Line(point1=(x_ct+x3[0] ,y_ct+y3[0] ), point2=(x_ct+x2[1],y_ct+y2[1])) 26 s1.Line(point1=(x_ct+x2[1] ,y_ct+y2[1] ), point2=(x_ct+x3[1],y_ct+y3[1])) 27 s1.Line(point1=(x_ct+x3[1] ,y_ct+y3[1] ), point2=(x_ct+x2[2],y_ct+y2[2])) 28 s1.Line(point1=(x_ct+x2[2] ,y_ct+y2[2] ), point2=(x_ct+x3[2],y_ct+y3[2])) 29 s1.Line(point1=(x_ct+x3[2] ,y_ct+y3[2] ), point2=(x_ct+x2[3],y_ct+y2[3])) 30 s1.Line(point1=(x_ct+x2[3] ,y_ct+y2[3] ), point2=(x_ct+x3[3],y_ct+y3[3])) 31 s1.Line(point1=(x_ct+x2[3] ,y_ct+y2[3] ), point2=(x_ct+x3[4],y_ct+y3[4])) 32 s1.Line(point1=(x_ct+x3[4] ,y_ct+y3[4] ), point2=(x_ct+x2[4],y_ct+y2[4])) 33 s1.Line(point1=(x_ct+x2[4] ,y_ct+y2[4] ), point2=(x_ct+x3[5],y_ct+y3[5])) 34 s1.Line(point1=(x_ct+x3[5] ,y_ct+y3[5] ), point2=(x_ct+x2[5],y_ct+y2[5])) 35 s1.Line(point1=(x_ct+x2[5] ,y_ct+y2[5] ), point2=(x_ct+x3[6],y_ct+ y3[6])) 36 p = mdb.models[ Model−1 ]. parts [ Part−1 ] 37 f = p. faces 38 pickedFaces = f.getSequenceFromMask(mask=(’[#a ]’, ), ) 39 e, d1 = p.edges, p.datums 40 p. PartitionFaceBySketch ( faces=pickedFaces, sketch=s) 41 s.unsetPrimaryObject () 42 del mdb.models[ Model−1 ].sketches [ \_\ _profile \_\_ ] Script 3.8 Step 05 in Python. 34 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements • Step 06: Mesh seeds. Mesh seeds are placed on the crown segments in order to make each sub-part coincide with a triangular element. To do so, each segment must coincide with an edge, without the possibility to let its number change. Hence, we will meshed with ’tri free’conditions. Mesh module → Seed Edges → Choose the crown edges → Done → Basic: (By number: Number of Elements: 1; Bians: None); Constraints: ’Do not allow the number of elements to change’. Figure 3.20 Mesh seeds of the crowns segments 06. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2e = p.edges 3pickedEdges = e.getSequenceFromMask(mask=(’[#e1ffffff #ffff ]’, ), ) 4p.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FIXED) Script 3.9 Step 06 in Python. • Step 07: Mesh seeds. Mesh seeds for b cases are placed on the oblique edges of the half-hexagon, with the values explained in the previous section. Mesh → Seed Edges → Choose the oblique halfhexagon edges → Done → Basic: (By number: Number of Elements: number of elements on the outer circumference; Bias: None); Constraints: Állow the number of elements to increase or decrease’. Figure 3.21 Mesh seeds of the half-hexagon oblique edges. 3.3 Mesh generation 35 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2e = p.edges 3pickedEdges = e.getSequenceFromMask(mask=(’[#8000000 #20000 ]’, ), ) 4p.seedEdgeByNumber(edges=pickedEdges, number= number of elements in the last crown) Script 3.10 Step 07 in Python. • Step 08: Mesh seeds. Remaining mesh seeds for b cases are placed on the horizontal edges of the half-hexagon and the vertical edge inside ofit, with the values explained in the previous section. Mesh → Seed Edges → Choose the two horizontal hexagon edges and the three segments between the half hexagon and the outer crown → Done → Basic: (By number: Number of Elements: ceiling of the result of the division of the number of elements on the outer circumference divided by three; Bias: None); Constraints: “Allow the number of elements to increase or decrease”. Figure 3.22 Mesh seeds for the remaining segments.. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2e = p.edges 3pickedEdges = e.getSequenceFromMask(mask=(’[16000000 50000 ]’, ), ) 4p.seedEdgeByNumber(edges=pickedEdges, number= ceil(number of elements in the last crown / 3)) Script 3.11 Step 08 in Python. • Step 09: Mesh seeds. Mesh seeds. Let us impose an approximate velue of global seeds in order not to let the size of the crown elements determine the average size of the elements. Out of the hexagon and far away from it, everything should remain approximately equal for every study. Mesh module → Seed Part → (Approximate global size: 0.05, tick on curvature control: Maximum deviation factor: 0.1, Minimum size control: by fraction of global size: 0.1). 36 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements Figure 3.23 Global seeds. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2p. seedPart ( size =0.05, deviationFactor =0.1, minSizeFactor=0.1) Script 3.12 Step 09 in Python. • Step 10: Mesh control. After several tries, a free triangular mesh confers symmetry to the model outside the half-hexagon. The structured meshes also producesasymmetric outcome, but triangles become further from equilateral as we distance from the half-hexagon. That is not a positive trait for a mesh. Mesh Module → Assign Mesh Controls → Select the whole part → (Element Shape: Tri, Technique: Free). Figure 3.24 Step 8. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2f = p. faces 3.3 Mesh generation 37 3pickedRegions = f.getSequenceFromMask(mask=(’[#3fffffff ]’, ), ) 4p.setMeshControls( regions=pickedRegions, elemShape=TRI) Script 3.13 Step 10 in Python. • Step 11: Part meshing..Mesh Module → Mesh Part → Done. Between the half-hexagon and the crowns, we obtain eitherasymmetric mesh or a mesh in which both halves differ in one vertex connection, that is to say, the orientation of two elementsédges. Inside the crowns, each sub-part is a different element, hence, we always obtain the desired mesh. Figure 3.25 Step 8. Figure 3.26 Step 8. 1p = mdb.models[ Model−1 ]. parts [ Part−1 ] 2p.generateMesh() Script 3.14 Step 11 in Python. • Step 12: Naming, saving and exporting.. We turn the values of the parameters into strings to compose the naming code. Then, the changing values are introduced in their corresponding spot in the code. We pick only the first 4 or 5 digits of each α value for the code, with the first digit always being a 0 (as a mantissa). 38 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 1a = mdb.models[ Model−1 ].rootAssembly 2a.DatumCsysByDefault(CARTESIAN) 3p = mdb.models[ Model−1 ]. parts [ Part−1 ] 4a. Instance (name= Part−1−1 , part=p, dependent=ON) 5alpha1_s=str(alpha1).replace(’.’, ”) 6alpha2_s=str(alpha2).replace(’.’, ”) 7alpha3_s=str(alpha3).replace(’.’, ”) 8tipo1_s=str(tipo1)) 9tipo2_s=str(tipo2)) Script 3.15 Step 12 in Python. 1!\ ctext [RGB]{252,190,17}{nofelements\_s=str(nofelements))}! 2nombre=nofelements\_s+"E\_tipo\_"+tipo1\_s+"\_tipo \_"+tipo2\_s+"\_"+alpha1\_s[0:5]+"\_"+alpha2\_s[0:4]+"\_"+alpha3\ _s[0:4] 3mdb.Job(name=nombre, model= Model−1 , description= , type=ANALYSIS, 4atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 5memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 6explicitPrecision =DOUBLE_PLUS_PACK, nodalOutputPrecision=FULL, 7echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint =OFF, 8userSubroutine= , scratch = , resultsFormat =ODB) 9mdb.jobs[nombre].writeInput (consistencyChecking=OFF) 10 mdb.saveAs( 11 pathName= C:/Path/ +nofelements_s+ E/ +nombre) Script 3.16 Step 13 in Python. 3.4 Type5crowns A minor change must be introduced for type 5 crowns, which cannot be built fromacircumference due to the vertical line. In this case, the crown should be cut by using the following segment arrangement: 1s.ArcByCenterEnds(center=(x_ct,y_ct), point1=(x_ct−rco2, y_ct+0), point2=(x_ct+x3[nofelements+1], y_ct+y3[ nofelements+1]), direction =CLOCKWISE) 2s.ArcByCenterEnds(center=(x_ct,y_ct), point1=(x_ct+rco2, y_ct) , 3point2=(x_ct+x3[nofelements−1],y_ct+y3[nofelements−1]), direction =COUNTERCLOCKWISE) 4s.Line(point1=(x_ct+x3[nofelements+1], y_ct+y3[nofelements+1]), point2=(x_ct+x3[nofelements−1], y_ct+y3[nofelements −1])) Script 3.17 Type 5 crown definition in Python. In order to achieve the desired result: izd dch vf vi vc vacvpc cnt Figure 3.27 Type 5 crown definition. 3.5 Process Improvement Making a step-by-step process like this required identifying errors separately at each line. That is why this method was chosen. However, after having more experience debugging the code, the process can be time-optimized by doing as follows: instead ofinterchanging different steps that draw segments and erase 3.6 Parametric sweep 39 regions, using only one step to make all cuts and then remove the SER. The result would make all the macros more interchangeable and save some time. The meshing steps and the result would not change. 3.6 Parametric sweep This study analyzes: • The effect of the number of elements around the crack tip. • Effect ofacrown. • Effect ofasecond crown. • The effect of the size proportion between the singular elements and the crowns elements. • The effect of the mesh type of crowns. The following diagram describes the parameter sweep. 1Hexagon_sidelength = 0.1; % 10% of the shell s sidelength 2for Number_of_Elements=[4, 6, 8, 12, 24] 3for crown_type_1 = [1, 2, 3, 4, 5]: 4for crown_type_2 = [0, 1, 2, 3, 4, 5]; 5for alpha1 =0.05:0.05: max_allowed alpha1_size 6for alpha2=min_allowed_alpha2_size_0.05:max_allowed alpha2_size 7for alpha3=min_allowed_alpha3_size_0.05:max_allowed alpha3_size 8Create part 9Create half hexagon % In this study, it is always the same length , it depends on the parameter hexagon_sidelength . 10 Coordinates (connected by lines ): % A unitary hexagon s vertices coordinates times side−length 11 x_hexagon = x_crack_tip +[1,0.5,0,−0.5,−1]* hexagon_sidelength 12 y_hexagon = y_crack_tip +[0, sqrt (3) /2, sqrt (3) /2, sqrt (3) /2,0]* hexagon_sidelength 13 Create auxiliary lines by segment coordinates : 14 Coordinates (connected by lines ) 15 x_aux = x_crack_tip +[0.5,1* hexagon_sidelength ]; x_crack_tip+[−0.5,−1*hexagon_sidelength ]; x_crack_tip +[0, 0] 16 y_aux = y_crack_tip +[1, sqrt (3)/2*hexagon_sidelength ]; y_crack_tip +[1, sqrt (3)/2*hexagon_sidelength ]; y_crack_tip +[1,0] 17 Create circumferences 18 Coordinates (by center and radius ): % Only the first two ones if there is only one crown 19 inner : 20 x_center = x_crack_tip , y_center = y_crack_tip ; radius = alpha1*hexagon_sidelength 21 medium: 22 x_center = x_crack_tip , y_center = y_crack_tip ; radius = alpha2*hexagon_sidelength 23 outer : 24 x_center = x_crack_tip , y_center = y_crack_tip ; radius = alpha3*hexagon_sidelength 25 Hollow inner circumference 26 Create crown lines 27 Coordinates: discussed in annex % that depend on alpha_1, alpha_2, alpha_3 and type_1 and type_2 28 Mesh seeds 29 Mesh control 30 Mesh part 31 Save ". inp" 32 end 33 end 34 end 35 end 36 end 37 end Script 3.18 MATLAB pseudo-code. The process could theoretically employ always the same commands withafew lines changed ifit was not for the internal numeration of the Abaqus elements (edges, vertices, regions, etc.), which differs for every different cut. For instance, the identification code assigned to the singular elements sub-part does not coincide in each case study; so when the region is erased, we will occasionally delete other sub-parts (the crown or exterior regions). However, for the first crown types 1, 2 and 3, and second crown types 1, 2, and 3; we can (most times) interchange sections of the code. 46 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 3.9 Simulation For the case simulation we will be usingaslightly modified version of [3], [5], [17]. Each of the program’s sections is discussed in its comments. The new program additions are the automatic imposition of boundary conditions, which is discussed in depth. Without it, the analysis of such an amount of meshes would not be feasible. 1%% Main loop. It must be in the same folder as the folders containing the . inp files . 2 3clear all ,close all ,clc 4 5%% This program solves our case study and gives displacements in the nodes by FEM method as outcome. 6 7%% Identification of the folder with . inp files . 8id= 4E_tipo_1_tipo_0 ;%% Name of the folder 9 10 %% Input (write values) : 11 hex_sidelength =0.1; % Half−hexagon side−length. 12 typeelement=1; % 1: with singular elements // 0: without singular elements 13 vct=[0.00 0.00 0.00]; % Vector with crack tip coordinates . 14 G=1; k_robin=1; NUM=400; 15 16 %% Initialize variables . 17 [name_index, file_index ,alpha1_index,alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , singular_FE_sidelength_index ,new_entrance_inp_index] = step_01_set_variables () ; 18 19 %% Obtener parámetros del nombre 20 [nofelements , nofcrowns, tipo1 , tipo2 , separadores ] = step_02_obtain_parameters (id) ; 21 22 %% Name index 23 [alpha1_index,alpha2_index,alpha3_index, singular_FE_sidelength_index , file_index ,name_index,ncasos] = step_03_name_archives(nofelements ,nofcrowns,tipo1 , tipo2 , hex_sidelength , id ,name_index, file_index , singular_FE_sidelength_index ,alpha1_index,alpha2_index,alpha3_index); 24 25 %% Loop to analyze every case study 26 for s=1:size (name_index,1) 27 close all ; 28 29 %% Change in the work directory (change it to fit the PC in use) 30 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 31 32 %% Specifying the values of certain parameters to fit the case index 33 alpha1=alpha1_index (:, s); 34 SE_sidelength=singular_FE_sidelength_index (s); 35 file = file_index (s ,:) ; 36 inp_malla= fileread ( file ); 37 38 %% Obtain meshes 39 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\Codigo Matlab ]) ; 40 [coordenadas,coordenadasN5,coordenadas5, triangulos , triangulos5 ]=step_04_mesher(inp_malla,nofelements,SE_sidelength, alpha1,typeelement, vct) ; 41 42 %% Specifying boundary conditions 43 [neumann,neumann5,dirichlet, dirichlet5 , robin , robin5,nodos_calculo]=step_05_boundary_conditions(coordenadas, vct , triangulos5 ) ; 44 45 %% Counting elements 46 [numNodos,numTriangulos,numTriangulos5,numLadosNeu,numLadosNeu5,numLadosDir,numLadosDir5,numLadosRob, numLadosRob5]=step_06_count_elements(coordenadas,triangulos,triangulos5,neumann,neumann5,dirichlet,dirichlet5, robin,robin5); 47 48 %% K_mn matrix of the singular elements 49 RobMat=step_07_robin_matrix(); 50 51 %% Variable initialization 52 [K,b,KROBIN,u] = step_08_initialize_system_matrixes (numNodos); 53 54 %% Stiffness matrix of the system 55 K = step_09_stiffness_matrix (coordenadas, triangulos , triangulos5 ,numTriangulos,numTriangulos5,G,K); 3.9 Simulation 47 56 57 %% Generation of matrices to assamble alongside stiffness matrix due to Robin conditions 58 KROBIN = step_10_robin_stiffness_matrix(numLadosRob,numLadosRob5,robin,robin5,coordenadas,k_robin,RobMat, numTriangulos,triangulos,KROBIN,NUM); 59 60 %% Assemble Robin conditions in the stiffness matrix: 61 K=K+KROBIN; 62 63 %% Neumann conditions 64 b = step_11_neumann_conditions(b,numLadosNeu,numLadosNeu5,neumann,neumann5,coordenadas,numTriangulos,RobMat, numTriangulos5,triangulos,triangulos5); 65 66 %% Dirichlet conditions 67 [b,Dir,Dir5] = step_12_dirichlet_conditions (b,numLadosDir,numLadosDir5,dirichlet, dirichlet5 ,coordenadas, numTriangulos,RobMat,numTriangulos5,triangulos, triangulos5 ); 68 69 %% Independent vector 70 b=b−K*u; 71 72 %% Solution 73 nodosDirTotales=[Dir;Dir5]; %forma una matriz con todos los nodos con condiciones Dirichlet 74 nodosLibres= setdiff (1:numNodos,nodosDirTotales); %excluye a los nodos con condiciones Dirichlet 75 u(nodosLibres)=K(nodosLibres,nodosLibres)\b(nodosLibres); 76 77 %% Displacement along the crack 78 desplazamientos_calculo =[ desplazamientos_calculo ; u(nodos_calculo) ]; 79 coordenadas_calculo=[coordenadas_calculo ; coordenadas(nodos_calculo ,1) ]; 80 end 81 82 %% Regression 83 [ ajuste , ajuste_2 , slope_index ,max_2nd_crown_change,casos_buenos] = step_13_numerical_regressions(alpha1_index, alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , tipo1,tipo2 , id, separadores ) ; 84 85 %% Saving 86 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 87 save(id); 88 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\Codigo Matlab ]) ; 89 90 end Script 3.22 Step 00 in MATLAB, main body. 3.9.1 step_01_set_variables 1function [name_index, file_index ,alpha1_index,alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , singular_FE_sidelength_index ,new_entrance_inp_index] = step_01_set_variables () 2%% This code initializes vectors that will be used in the program. 3name_index=[]; 4file_index =[]; 5alpha1_index =[]; 6alpha2_index =[]; 7alpha3_index =[]; 8desplazamientos_calculo =[]; 9coordenadas_calculo =[]; 10 singular_FE_sidelength_index =[]; 11 new_entrance_inp_index=[]; 12 end Script 3.23 Step 01 in MATLAB. 3.9.2 step_02_obtain_parameters 1function [nofelements , nofcrowns, tipo1 , tipo2 , separadores ] = step_02_obtain_parameters (id) 2%% This code finds the crown type data by the naming code 3separadores= strfind (id , _); 48 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 4nofelements=str2num(id(1: separadores (1)−2)); 5if str2num(id( separadores (end)+1:end))~=0 6nofcrowns=2; 7else 8nofcrowns=1; 9end 10 tipo1=str2num(id( separadores (2)+1:separadores (3)−1)); 11 tipo2=str2num(id( separadores (4)+1:end)); Script 3.24 Step 02 in MATLAB. 3.9.3 step_03_name_archives A reduced version of the code is here presented. The complete code in in the appendix. The scripts is explained in the commented lines. 1 2function [alpha1_index,alpha2_index,alpha3_index, singular_FE_sidelength_index , file_index ,name_index,ncasos] = step_03_name_archives(nofelements ,nofcrowns,tipo1 , tipo2 , hex_sidelength , id ,name_index, file_index , singular_FE_sidelength_index ,alpha1_index,alpha2_index,alpha3_index, separadores ) 3cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 4new_entrance=[]; 5 6%% This code generates the name index that we will use to identify files and read them. It also generates an index of alpha_1, alpha_2 and alpha_3 by repeating Abaqus loops with the same values . Since their names were originally created in Abaqus, the strings generated in MATLAB will probably differ in the number of zeros added or in one decimal due to rounding. That is why some lines are added to fit the format to Abaqus . The limits are those indicated previously for each crown type. The process is as follows : we turn the alpha coefficients into strings with the most likely format. If that file is not found, we know that as much we have to add or subtract 0.0001, which is the format limit (four decimals precision ) . The configurations are several : alpha_2 can augment, decrease or remain equal, but if there is a second crown, alpha_3 can simultaneously do the same. All the possibilities are explored , and the first one to be found is chosen. 7 8%% 1 Crown: 9 10 %% type_1_type_0, type_2_type_0, type_3_type_0 11 if nofcrowns==1 && (tipo1==1 || tipo1==2 || tipo1==3) 12 alpha1_in=0.05; 13 alpha1_in_plus=0; 14 alpha2_in_n=1/cos(pi/nofelements); 15 alpha2_in_plus=0.05; 16 end 17 18 %% type_4_type_0 19 if nofcrowns==1 && tipo1==4 20 alpha1_in=0.05; 21 alpha1_in_plus=0; 22 alpha2_in_n=1/cos(pi /2/ nofelements); 23 alpha2_in_plus=0; 24 end 25 %% type_5_type_0 for 4, 6 and 8 elements. 26 if nofcrowns==1 && tipo1==5 && (nofelements==4 || nofelements==6 || nofelements==8) 27 alpha1_in=0.05; 28 alpha1_in_plus=0; 29 alpha2_in_n=1/cos(pi /2/ nofelements); 30 alpha2_in_plus=0.0001; 31 end 32 33 %% Two crowns 34 35 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && (tipo2==1 || tipo2==2 || tipo2==3) 36 alpha1_in=0.05; 37 alpha1_in_plus=0; 38 alpha2_in_n=1/(cos(pi/nofelements)); 39 alpha2_in_plus=0.05; 40 alpha3_in_n=1/(cos(pi/nofelements)); 41 alpha3_in_plus=0.05; 42 end 43 3.9 Simulation 49 44 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && tipo2==4 45 alpha1_in=0.05; 46 alpha1_in_plus=0; 47 alpha2_in_n=1/(cos(pi/nofelements)); 48 alpha2_in_plus=0; 49 alpha3_in_n=1/(cos(pi /2/ nofelements)); 50 alpha3_in_plus=0; 51 end 52 53 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && tipo2==5 54 alpha1_in=0.05; 55 alpha1_in_plus=0; 56 alpha2_in_n=1/(cos(pi/nofelements)); 57 alpha2_in_plus=0.05; 58 alpha3_in_n=1/(cos(pi /2/ nofelements)); 59 alpha3_in_plus=0.05; 60 end 61 62 %% One−crown 63 if nofcrowns==1, 64 for alpha1=alpha1_in+alpha1_in_plus :0.05:0.85 65 for alpha2=alpha1∗alpha2_in_n+alpha2_in_plus :0.05:0.85 66 singular_FE_sidelength =alpha1∗hex_sidelength; 67 alpha1_s=num2str((alpha1) ,4) ; 68 alpha1_s=strrep (alpha1_s, '.',' ' ) ; 69 alpha2_s=num2str((alpha2), '%2.15f'); 70 alpha2_s=strrep (alpha2_s, '.',' ' ) ; 71 alpha2_s=alpha2_s(1:4) ; 72 Files=dir ([ '∗' ,'_',alpha1_s, '_',alpha2_s(1:end), ' ∗.inp ']) ; 73 if isempty( Files ) 74 alpha2_s=strrep (num2str(alpha2+0.0001,'%2.15f') , '.',' ' ); 75 alpha2_s=alpha2_s(1:4) ; 76 Files =[]; 77 Files=dir ([ '∗' ,'_',alpha1_s, '_',alpha2_s(1:end), ' ∗.inp ']) ; 78 if isempty( Files ) 79 alpha2_s=strrep (num2str(alpha2−0.0001,'%2.15f') , '.',' ' ); 80 alpha2_s=alpha2_s(1:4) ; 81 end 82 end 83 new_entrance=[id alpha1_s '_'alpha2_s ]; 84 new_entrance_inp=[id '_'alpha1_s '_'alpha2_s '.inp ']; 85 name_index=strvcat(name_index, new_entrance); 86 file_index = strvcat ( file_index ,new_entrance_inp); 87 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; singular_FE_sidelength ]; 88 alpha1_index=[alpha1_index,alpha1 ]; 89 alpha2_index=[alpha2_index,alpha2 ]; 90 end 91 end 92 end 93 94 %% Two−crowns 95 if nofcrowns==2 96 for alpha1=alpha1_in+alpha1_in_plus :0.05:0.85 97 for alpha2=alpha1∗alpha2_in_n+alpha2_in_plus :0.05:0.85 98 for alpha3=alpha2∗alpha3_in_n+alpha3_in_plus :0.05:0.85 99 % Singular element side−length 100 singular_FE_sidelength =alpha1∗hex_sidelength; 101 102 % alpha1 103 alpha1_s=num2str((alpha1) ,4) ; 104 alpha1_s=strrep (alpha1_s, '.',' ' ) ; 105 106 % alpha2 107 alpha2_s=num2str((alpha2), '%2.15f'); 108 alpha2_s=strrep (alpha2_s, '.',' ' ) ; 109 alpha2_s=alpha2_s(1:4) ; 110 111 % alpha3 112 alpha3_s=num2str((alpha3), '%2.15f'); 113 alpha3_s=strrep (alpha3_s, '.',' ' ) ; 114 if length(alpha3_s)>3, alpha3_s=alpha3_s(1:4) ; else , alpha3_s=alpha3_s(1:end); end 50 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 115 116 % = alpha2, = alpha3 117 Files=dir ([ *,_,alpha1_s, _,alpha2_s(1:end), _,alpha3_s(1:end), *.inp ]) ; 118 if isempty( Files ), 119 % = alpha2, + alpha3. 120 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 121 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 122 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 123 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 124 alpha3_saux=alpha3_saux(1:4); 125 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 126 if ~isempty( Files ) , 127 alpha2_s=alpha2_saux; 128 alpha3_s=alpha3_saux; 129 else 130 % = alpha2, −alpha3. 131 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 132 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 133 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 134 alpha3_saux=alpha3_saux(1:4); 135 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 136 if ~isempty( Files ) , 137 alpha2_s=alpha2_saux; 138 alpha3_s=alpha3_saux; 139 else 140 % + alpha2, = alpha3 141 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 142 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 143 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 144 alpha2_saux=alpha2_saux(1:4); 145 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 146 if ~isempty( Files ) , 147 alpha2_s=alpha2_saux; 148 alpha3_s=alpha3_saux; 149 else 150 % + alpha2, + alpha3 151 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 152 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 153 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 154 alpha2_saux=alpha2_saux(1:4); 155 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 156 alpha3_saux=alpha3_saux(1:4); 157 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 158 if ~isempty( Files ) , 159 alpha2_s=alpha2_saux; 160 alpha3_s=alpha3_saux; 161 else 162 % + alpha2, − alpha3 163 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 164 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 165 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 166 alpha2_saux=alpha2_saux(1:4); 167 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 168 alpha3_saux=alpha3_saux(1:4); 169 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 170 if ~isempty( Files ) , 171 alpha2_s=alpha2_saux; 172 alpha3_s=alpha3_saux; 173 else 174 % − alpha2, = alpha3 175 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 176 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 177 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 178 alpha2_saux=alpha2_saux(1:4); 179 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 180 if ~isempty( Files ) , 181 alpha2_s=alpha2_saux; 182 alpha3_s=alpha3_saux; 183 else 184 % − alpha2, + alpha3 185 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 3.9 Simulation 51 186 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 187 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 188 alpha2_saux=alpha2_saux(1:4); 189 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 190 alpha3_saux=alpha3_saux(1:4); 191 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 192 if ~isempty( Files ) , 193 alpha2_s=alpha2_saux; 194 alpha3_s=alpha3_saux; 195 else 196 % − alpha2, − alpha3 197 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 198 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 199 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 200 alpha2_saux=alpha2_saux(1:4); 201 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 202 alpha3_saux=alpha3_saux(1:4); 203 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 204 if ~isempty( Files ) , 205 alpha2_s=alpha2_saux; 206 alpha3_s=alpha3_saux; 207 end 208 end 209 end 210 end 211 end 212 end 213 end 214 end 215 end 216 if isempty( Files ) 217 new_entrance_inp=[]; 218 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; []]; 219 name_index=strvcat(name_index, new_entrance); 220 alpha1_index=[alpha1_index ,[]]; 221 alpha2_index=[alpha2_index ,[]]; 222 alpha3_index=[alpha3_index ,[]]; 223 else 224 new_entrance_inp=[id _alpha1_s _alpha2_s _alpha3_s .inp ]; 225 name_index=strvcat(name_index, new_entrance); 226 file_index = strvcat ( file_index ,new_entrance_inp); 227 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; singular_FE_sidelength ]; 228 alpha1_index=[alpha1_index,alpha1 ]; 229 alpha2_index=[alpha2_index,alpha2 ]; 230 alpha3_index=[alpha3_index,alpha3 ]; 231 end 232 end 233 end 234 end 235 end 236 237 ncasos=length( file_index ); Script 3.25 Step 03 in MATLAB. 3.9.4 step_04_mesher 1function [coordenadas,coordenadasN5,coordenadas5, triangulos , triangulos5 ] = step_04_mesher(inp_malla,nofelements, sidelength ,alpha,typeelement, vct); 2xcentro=vct(1) ; ycentro=vct(2) ; zcentro=vct(3) ; 3 4%% Coordinates matrix 5begin_reading_coordenadas= strfind (inp_malla, *Node ); % Reads from "*Node". 6end_reading_coordenadas= strfind (inp_malla, *Element ) ; % to "*Element". 7coordenadas_aux=[inp_malla(begin_reading_coordenadas+8:end_reading_coordenadas−2)]; % The . inp coordinate file is reduce to the necessary part , the one corresponding to the nodal coordinates . The +5 is due to read the beginning of the numbers, and not the heading. The −2 is bound to not read the following heading characters . 52 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 8coordenadas_aux= strsplit (coordenadas_aux, \n ); % The text is divided in lines by the newline character . 9coordenadas_aux=coordenadas_aux ; % Vertically arranged text . 10 for i=1:length(coordenadas_aux) 11 coordenadas(i ,:) =str2num(cell2mat(coordenadas_aux(i))); % Every row is turned into a text array , and every text array (as a cell ) is turned into a numeric value. 12 end 13 coordenadas (:,1) =[]; % The identification number column is deleted . 14 15 %% Elements connectivity matrix 16 17 begin_reading_triangulos = strfind (inp_malla, *Element, ) ; % The beginning of the Triangle matrix is found 18 end_reading_triangulos = strfind (inp_malla, *End Part ) ; % The end of the Triangle matrix is found 19 triangulos_aux =[inp_malla( begin_reading_triangulos +21: end_reading_triangulos −2)]; % Their values are introduced in a variable 20 triangulos_aux = strsplit ( triangulos_aux , \n ); % The text is divided in lines by the newline character . 21 triangulos_aux =triangulos_aux ; % The matrix is trasposed to get the correct format. 22 for i=1:length( triangulos_aux ) 23 triangulos (i ,:) =str2num(cell2mat( triangulos_aux ( i))); % The triangle variable line is turned into numbers 24 end 25 triangulos (:,1) =[]; % The first column gets deleted since it is unnecessary 26 27 %% Singular−Elements Coordinates Matrix 28 29 N=2*nofelements; % The number of regular polygon edges whose upper half will be singular elements. 30 n=(0:nofelements); 31 r=sidelength *1/2; % alpha=1/2, the position of the node in the singular −behavioured edge. 32 if typeelement==1 % If there are singular elements. 33 x(n+1)=r*cos(2*pi*n/N); 34 y(n+1)=r*sin(2*pi*n/N); % Their x and y regular polygon coordinates are calculated 35 z(n+1)=zcentro; % Introduce z value. 36 coordenadas_singulares =[xcentro , ycentro , zcentro ; x , y , z ]; 37 elseif typeelement==0 % If there are no singular elements 38 coordenadas_singulares =[xcentro , ycentro , zcentro ]; 39 end 40 coordenadas_singulares (:,3) =[]; % The third column gets deleted , since it is unnecessary. 41 42 %% Singular elements connectivity matrix. 43 44 r=sidelength ; 45 x(n+1)=r*cos(2*pi*n/N); 46 y(n+1)=r*sin(2*pi*n/N); 47 z(n+1)=0; % Introduce z coordinate . 48 for n=1:nofelements+1 49 closestpoint (n)=dsearchn(coordenadas,[x(n) y(n)]) ; % This command finds the closest point (in the mesh) to the ones obtained by our formula (they may not perfectly coincide due to minimal Abaqus−MATLAB errors). 50 end 51 if typeelement==1 % The connectivity matrix is created for singular elements. The points are ordered clockwise: the first node is the central one, the second to its right , and the third one is the pre− existent one ( already in the mesh, since it is part of both the singular and non−singular element), the fourth point will be the other pre− existent node and the last one will be the remaining intermediate point . 52 for i=1:nofelements 53 triangulos5 (i ,:) =[size (coordenadas,1)+1 size (coordenadas,1)+1+i closestpoint (i) closestpoint (i+1) size ( coordenadas,1)+2+i]; 54 end 55 elseif typeelement==0 56 for i=1:nofelements % In this case, no halfway points exist . 57 triangulos5 (i ,:) =[size (coordenadas,1)+1 closestpoint (i) closestpoint (i+1)]; 58 end 59 triangulos =[ triangulos ; triangulos5 ]; 60 triangulos5 =[]; 61 end 62 %% This line assembles the coordinate matrix. 63 coordenadas5=coordenadas_singulares ; 64 coordenadasN5=coordenadas; 65 coordenadas=[coordenadas; coordenadas_singulares ]; Script 3.26 Step 04 in MATLAB. 3.9 Simulation 53 3.9.5 step_05_boundary_conditions Function step_05_boundary_conditions part 1 1%% This line finds the plate edges. 2[ Lista_nodos_borde_izquierdo , Lista_nodos_borde_superior , ... 3Lista_nodos_borde_derecho, Lista_nodos_borde_inferior , ... 4Nodos_borde_izquierdo, Nodos_borde_superior, ... 5Nodos_borde_derecho, Nodos_borde_inferior ], ... 6=encuentra_bordes(coordenadas, tolerancia ); Script 3.27 Step 05 in MATLAB. Function encuentra_bordes 1function [ Lista_nodos_borde_izquierdo , Lista_nodos_borde_superior ,Lista_nodos_borde_derecho, Lista_nodos_borde_inferior ,Nodos_borde_izquierdo, Nodos_borde_superior,Nodos_borde_derecho, Nodos_borde_inferior]=encuentra_bordes(coordenadas, tolerancia ) 2 3%% This function generates 8 lists : 4 with the nodal number of each edge nodes (upper, lower, right and left ) , 4 with the coordinates of said nodes. 4 5%% Initialize matrices 6Nodos_borde_derecho=[]; Lista_nodos_borde_derecho=[]; 7Nodos_borde_izquierdo=[]; Lista_nodos_borde_izquierdo =[]; 8Nodos_borde_superior=[]; Lista_nodos_borde_superior =[]; 9Nodos_borde_inferior =[]; Lista_nodos_borde_inferior =[]; 10 11 %% Obtain the boundary points from our node list 12 frontera =boundary(coordenadas (:,1) ,coordenadas (:,2) ); 13 14 %% Obtain the upper and lower (y coordinate ) and most to the right and to the left limits of the coordinates (edges ) 15 coordenada_borde_derecho=max(coordenadas(:,1)); 16 coordenada_borde_izquierdo=min(coordenadas (:,1) ); 17 coordenada_borde_superior=max(coordenadas(:,2)); 18 coordenada_borde_inferior =min(coordenadas (:,2) ); 19 20 %% We define a tolerance in case there is a mild discordance among the points coordinates 21 tol=1e−10; 22 23 %% Compute borders. We add to each list the nodes whose coordinates fit each limit . If a nodes coordinates coincide with the right or left limits ( including tolerance ) , they are added to the list . 24 for i=1:length( frontera ) 25 nodo=frontera(i) ; 26 if coordenadas(nodo,1)>=coordenada_borde_derecho−tol && coordenadas(nodo,1)<=coordenada_borde_derecho+tol, 27 Nodos_borde_derecho=[Nodos_borde_derecho; coordenadas(nodo,:)]; 28 Lista_nodos_borde_derecho=[Lista_nodos_borde_derecho; nodo]; 29 end 30 if coordenadas(nodo,1)>=coordenada_borde_izquierdo−tol && coordenadas(nodo,1)<=coordenada_borde_izquierdo+tol, 31 Nodos_borde_izquierdo=[Nodos_borde_izquierdo; coordenadas(nodo,:) ]; 32 Lista_nodos_borde_izquierdo=[Lista_nodos_borde_izquierdo ; nodo]; 33 end 34 if coordenadas(nodo,2)>=coordenada_borde_superior−tol && coordenadas(nodo,2)<=coordenada_borde_superior+tol, 35 Nodos_borde_superior=[Nodos_borde_superior; coordenadas(nodo ,:) ]; 36 Lista_nodos_borde_superior=[Lista_nodos_borde_superior ; nodo]; 37 end 38 if coordenadas(nodo,2)>=coordenada_borde_inferior−tol && coordenadas(nodo,2)<=coordenada_borde_inferior+tol, 39 Nodos_borde_inferior=[Nodos_borde_inferior ; coordenadas(nodo,:) ]; 40 Lista_nodos_borde_inferior =[ Lista_nodos_borde_inferior ; nodo]; 41 end 42 end 43 44 %% Deletes repeated values with a call to elimina_dobles . 45 [ Lista_nodos_borde_izquierdo ,Nodos_borde_izquierdo]=elimina_dobles( Lista_nodos_borde_izquierdo , Nodos_borde_izquierdo); 46 [ Lista_nodos_borde_superior ,Nodos_borde_superior]=elimina_dobles( Lista_nodos_borde_superior ,Nodos_borde_superior); 47 [Lista_nodos_borde_derecho,Nodos_borde_derecho]=elimina_dobles(Lista_nodos_borde_derecho,Nodos_borde_derecho); 54 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 48 [ Lista_nodos_borde_inferior , Nodos_borde_inferior]=elimina_dobles( Lista_nodos_borde_inferior , Nodos_borde_inferior); 49 end Script 3.28 Function Encuentra_bordes. Function elimina_dobles 1function [Lista_nueva ,nodos_nuevos] = elimina_dobles ( Lista ,nodos) 2 3%% This function deletes repeated adjacent values to make a simple list of boundary nodes. 4Lista_nueva=Lista ; 5nodos_nuevos=nodos; % List of nodes and coordinates copied in a new vector 6eliminar =[]; % The pointer vector that will give us the values to erase is set . 7for i =1:(( length (unique(Lista_nueva)))−1) % If a value coincides with the next value in the array , its position is introduced in the vector . 8if Lista_nueva(i)==Lista_nueva(i+1) 9eliminar=[eliminar , i ]; 10 end 11 end 12 if Lista_nueva(end)==Lista_nueva(1) % A coincidence between the first and last is also checked. 13 eliminar=[eliminar length (Lista_nueva) ]; 14 end 15 Lista_nueva( eliminar ) =[]; 16 nodos_nuevos(eliminar ,:) =[]; % The repeated values are deleted . 17 end Script 3.29 Function elimina dobles. Function step_05_boundary_conditions part 2 1%% This section orders the edges nodes clockwise from the lower left node. 2[Lista_nodos, Lista_nodos_borde_izquierdo ,... 3Lista_nodos_borde_superior ,Lista_nodos_borde_derecho ... 4...=]roirefni_edrob_sodon_atsiL 5ordena_bordes( Lista_nodos_borde_izquierdo , ... 6Lista_nodos_borde_superior ,Lista_nodos_borde_derecho, ... 7Lista_nodos_borde_inferior ,Nodos_borde_izquierdo, ... 8Nodos_borde_superior, Nodos_borde_derecho, ... 9Nodos_borde_inferior) ; Script 3.30 Step 05 in MATLAB part 2. Function ordena_bordes 1%% The nodes in the vectors must be ordered clockwise, starting with the bottom left node. 2 3function [Lista_nodos, Lista_nodos_borde_izquierdo ,... 4Lista_nodos_borde_superior ,Lista_nodos_borde_derecho ... 5Lista_nodos_borde_inferior ]=ordena_bordes( Lista_nodos_borde_izquierdo , Lista_nodos_borde_superior ,... 6Lista_nodos_borde_derecho, Lista_nodos_borde_inferior ,Nodos_borde_izquierdo,Nodos_borde_superior ,... 7Nodos_borde_derecho,Nodos_borde_inferior) 8 9%% The right Edge is ordered from lower to higher according to their y coordinate . The right −edge nodes are introduced in said order , and then is flipped . 10 [nodos_sorted,nodos_orden]=sort(Nodos_borde_derecho(:,2)); 11 Lista_nodos_borde_derecho=Lista_nodos_borde_derecho(nodos_orden,:) ; 12 Lista_nodos_borde_derecho=flip (Lista_nodos_borde_derecho); 13 %% The left Edge does not require the flip . 14 [nodos_sorted,nodos_orden]=sort(Nodos_borde_izquierdo (:,2) ); 15 Lista_nodos_borde_izquierdo=Lista_nodos_borde_izquierdo(nodos_orden,:) ; 16 %% Idem for the top Edge nodes from left to right . 17 [nodos_sorted,nodos_orden]=sort(Nodos_borde_superior (:,1) ); 18 Lista_nodos_borde_superior=Lista_nodos_borde_superior(nodos_orden,:) ; 19 %% Finally the lower Edge receives the same treatment , but is flipped afterwards . 3.9 Simulation 55 20 [nodos_sorted,nodos_orden]=sort(Nodos_borde_inferior (:,1) ); 21 Lista_nodos_borde_inferior = Lista_nodos_borde_inferior (nodos_orden,:) ; 22 Lista_nodos_borde_inferior = flip ( Lista_nodos_borde_inferior ); 23 24 %% This section generates a list with all the nodes ordered , and the double values are deleted . 25 Lista_nodos=[Lista_nodos_borde_izquierdo ; Lista_nodos_borde_superior ; ... 26 Lista_nodos_borde_derecho; Lista_nodos_borde_inferior ]; 27 Nodos=[Nodos_borde_izquierdo; Nodos_borde_superior; Nodos_borde_derecho; ... 28 Nodos_borde_inferior ]; 29 [Lista_nodos,Nodos]=elimina_dobles(Lista_nodos,Nodos); 30 end Script 3.31 Function Ordena Bordes. Function step_05_boundary_conditions part 3 1%% For different boundary conditions , fill their definition here. 2%% Values must be introduced always clockwise. In our case: there are no nodes with Dirichlet boundary condition . Neumann nodes are top edge (value 1), lateral and bottom edge left half (value 0). Robin nodes are bottom edge right half (parameter k). Except for the singular element nodes, which are bottom central and the following two sideways. 3nodos_dirichlet =[]; 4nodos_neumann=[Lista_nodos_borde_izquierdo;... 5Lista_nodos_borde_superior ;... 6Lista_nodos_borde_derecho ;... 7Lista_nodos_borde_inferior ( floor (length ( Lista_nodos_borde_inferior ) /2)+3:end) ;... 8Lista_nodos_borde_izquierdo (1) ]; 9 10 %% We select the left , upper, right edges, and from the bottom edge, from the rightmost to the midway but three node. The left edge first node is added again to close the circle (between the last and first nodes). 11 if rem(length( Lista_nodos_borde_inferior ) ,2)==0 12 % Si el número de nodos en el borde inferior fuera par, los nodos con condición Robin ir ían desde el primero hasta la mitad menos dos nodos. 13 %% If the number of nodes in the bottom edge is even, Robin nodes span from the the first to the midway but two nodes. 14 nodos_robin=[ Lista_nodos_borde_inferior (1:( length ( Lista_nodos_borde_inferior ) /2)−2)]; 15 else 16 %% If the number is odd, the formula is slightly different : the downward nearest whole number, and subtract one to arrive to the desired node. 17 nodos_robin=[ Lista_nodos_borde_inferior (1: floor (length ( Lista_nodos_borde_inferior ) /2)−1)]; 18 end 19 20 %% No singular nodes affected by Dirichlet condition . 21 nodos_dirichlet5 =[]; 22 %% We apply the same reasoning to non−singular elements. In this case, Neumann nodes span from midway edge node to the following two left nodes. 23 if rem(length( Lista_nodos_borde_inferior ) ,2)==0 24 nodos_neumann5=[Lista_nodos_borde_inferior((length ( Lista_nodos_borde_inferior ) /2) :( length ( Lista_nodos_borde_inferior ) /2)+2) ]; 25 else 26 %% And again, we must differentiate if the number of elements in the bottom edge is even or odd. 27 nodos_neumann5=[Lista_nodos_borde_inferior(floor (length ( Lista_nodos_borde_inferior ) /2) : floor (length( Lista_nodos_borde_inferior ) /2)+2) ]; 28 end 29 nodos_robin5=[ Lista_nodos_borde_inferior ( floor (length ( Lista_nodos_borde_inferior ) /2) :−1: floor (length ( Lista_nodos_borde_inferior ) /2)−2) ]; 30 31 %% The nodes are again rearranged according to the nodal list in case they had changed during the process . In general , they may not maintain its order . 32 if ~isempty( nodos_dirichlet ) 33 nodos_dirichlet =reordena( nodos_dirichlet ,Lista_nodos); 34 else dirichlet =[]; 35 end 36 if ~isempty(nodos_neumann) 37 nodos_neumann=reordena(nodos_neumann,Lista_nodos); 38 else neumann=[]; 39 end 40 if ~isempty(nodos_robin) 41 nodos_robin=reordena(nodos_robin,Lista_nodos); 62 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 196 %% Curve plot 197 xv = linspace (min(x), max(x), 20); % Linear alpha_2 interpolation values 198 yv = linspace (min(y), max(y), 20); % Linear alpha_3 interpolation values 199 [X,Y] = meshgrid(xv, yv); % X−Y reticule generation 200 Z = griddata (x,y,z,X,Y); % Sparce data interpolation to fit the Surface. 201 surf (X,Y,Z) 202 end 203 204 %% Ajuste de los planos 205 Xcolv = x_pr (:) ; % alpha_2 from the level curve 206 Ycolv = y_pr (:) ; % alpha_3 from the level curve 207 Zcolv = z_pr (:) ; % alpha_1 from the alpha_1 index 208 Const = ones( size (Xcolv)); 209 Coefficients = [Xcolv Ycolv Const]\Zcolv; % This section computes the plane coefficients by the Least Squares Method. 210 XCoeff = Coefficients (1) ; 211 YCoeff = Coefficients (2) ; 212 CCoeff = Coefficients (3) ; 213 214 %% alpha_1−alpha_2−alpha_3 adjustment. 215 % Third figure : fitting planes (alpha_1 vs alpha_3 vs alpha_3) 216 figure (2) 217 L=plot3(x_pr,y_pr,z_pr, ro ); % Plot the original data points 218 set (L, Markersize ,2*get(L, Markersize )) % Making the circle markers larger 219 set (L, Markerfacecolor ,r)% Filling in the markers 220 hold on 221 [xx, yy]=meshgrid (0:0.1:1,0:0.1:1) ; % Generating a regular grid for plotting 222 zz = XCoeff * xx + YCoeff * yy + CCoeff; 223 surf (xx,yy,zz) % Plotting the surface 224 title (sprintf (z=(%f)*x+(%f)*y+(%f) ,XCoeff, YCoeff, CCoeff)) 225 xlabel (\alpha_2 ); ylabel (\alpha_3 ); zlabel (\alpha_1 ); 226 print (−dpng ,−r300 , [ graph_plano_ajuste_ , id , _.png ]) 227 ajuste_2=[−XCoeff −YCoeff 1 CCoeff]; %% Z coordinate pointing upwards 228 229 %% MATLAB Surface fit 230 Coeff_vectora=[0 0]; Rsquare_plane_a=0; 231 Coeff_vectorb=[0 0 0]; Rsquare_plane_b=0; 232 if size (x_pr,1)>=3 233 [ Coeff_vectora ,gofa]= fit ([x_pr, y_pr], z_pr, poly11 ); 234 Rsquare_plane_a=gofa.rsquare ; % Goodness of fit for the linear regression 235 end 236 if size (x_pr,1)>6 237 [Coeff_vectorb ,gofb]= fit ([x_pr, y_pr], z_pr, poly22 ); 238 Rsquare_plane_b=gofb.rsquare; % Goodness of fit for the parabolic regression 239 end 240 241 %% Right−solution−plane graphic 242 % First and second figures : surfaces (u_ct vs alpha_2 vs alpha_3 for each alpha_1) 243 figure (1) ; grid on; xlabel (\alpha_2 ); zlabel (u_1 ) ; ylabel (\alpha_3 ); hold on 244 [X,Y] = meshgrid (0:0.7:0.7,0.1:0.8:0.9) ; 245 Z = 2.25596*ones(size(X)); 246 surf (X,Y,Z); view(45,20); axis ([0 0.8 0 1 1.75 2.5]) ; 247 print (−dpng ,−r300 , [ graph_curva_ajuste_ , id , .png ]) ; 248 view([0 0 1]) ; print (−dpng ,−r300 , [ graph_curva_ajuste_arriba_ , id , .png ]) 249 250 %% Score criteria 251 if isempty(casos_buenos)==0 252 criterio_1 = length(casos_buenos)/ length(alpha1_index(change)); 253 theta = atan(mean(slope_index)); 254 criterio_2 = (pi/2+ theta ) /( pi /2) ; 255 criterio_3a = Rsquare_plane_a; 256 criterio_3b = Rsquare_plane_b; 257 criterio_4 =1−size(coordenadas,1) /1131; 258 nota_a=( criterio_1 + criterio_2 + criterio_3a + criterio_4 ) /4; 259 nota_b=( criterio_1 + criterio_2 + criterio_3b + criterio_4 ) /4; 260 end 261 if isempty(casos_buenos)==1 262 criterio_1 =0; criterio_2 =0; criterio_3a =0; criterio_3b =0; criterio_4 =0; 263 nota_a=0; nota_b=0; 264 end 265 end 3.10 Results 63 266 [nota_a nota_b]; 267 close all 268 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros b\Codigo Matlab ]) ; Script 3.37 Step 13 in MATLAB. 3.10 Results Once the meshes are obtained and simulated, the results must be analyzed in some manner. First, let us plot the results to haveaquick overview. One-crown cases: In each graphic, the number of singular elements and the crown type remains constant, there is one curve for every α1value, which relates the crack tip displacement (uct ) with α2. Figure 3.32 One possible typology: few intersections. Cut-off points are extrapolations (not included in the domain of the curve).. Figure 3.33 Another possible typology: several intersections.. Cut-off points of each curve with the horizontal segment (which is the exact solution) will give us the optimal proportion between radii alpha1 and alpha2 . And indirectly, the optimal proportion between singular elements and crown elements. So when the crown element size is plot versus the singular element given by the real solution, we obtain points that can be fit to a curve. Figure 3.34 With few instances, a fitting curve givesagood approximation. Figure 3.35 With more instances, the best-fitting curve resembles a parabola. Both curves will be used: the line for simplicity and the parabola for accuracy.. 64 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements Two-crown cases: in each graphic the number of singular elements and the crown types remains constant, and for each α1 there will beadifferent surface, which relates uct to α2 in one axis and α3 in the other axis. Since linear and parabolic adjustments are so close to each other, we will present parabolic adjustment graphics from now on. Figure 3.36 Typology: few intersections.. Figure 3.37 Typology: several intersections.. The 3D surfaces resemble an extruded version of the 2D curves, that is, their projection into an α2−α1 plane follow the same shape as the 2D cases. That is to say, the projections follow approximately the shape of the one-crown case curves. That indicates that the result has little dependency on the second crown, since independently from its size, the obtained result does not change considerably. The second crown effect is onlyaslight warp that turns the uct -level curves into parabolas bended towards the decreasing α2 direction (as shown in the image). This variation due to the third crown is measured by the code. Figure 3.38 Their projection resembles the previous result. α2 α3 uct Figure 3.39 Parabolic behavior induced by the second crown (effect in green). Now, the points that result in the proper result have three coordinates: α1 , α2 , α3 . When scattered they appear to haveaplanar distribution: As we can see, the results have a good fit with a plane whereas the one-crown cases have a better fit with a parabola. By now, we can appreciate that the two-crowns results fit adequately enough inaplane whereas the one-crown cases strongly fit in a parabola. However, it is worth noting that the points are not perfectly aligned in a plane either: they haveaparabolic tendency although less accentuated than one-crown cases. Projecting the points onto the α2−α1plane allows us to spot said a parabolic-like tendency. We can use this code to find outafirst approximation to what the best proportion between the sizes of the singular elements rosette and the first crown. Only the linear (for one crown) and planar (for two crowns) will be shown, as they are more easily understandable. The parabolic and quadratic-surface adjustment (just as their optimal proportion) will be shown and analysed only for a set of cases (a set of good configurations) 3.10 Results 65 Figure 3.40 Typology: few intersections.. Figure 3.41 Typology: several intersections.. which will be discussed later. For one-crown cases, we have (α2−α1) = p1α1+p0 . Since p0 will usually haveanear-zero value, p1 will approximately be ratio of the crown width divided by the rosette radius. For two-crown cases, we have α2=p2α3+p1α1+p0 . As stated before, the dependency on α3 is low ( p2 is near zero), and like before, p0 is near zero. Hence, again p2 is a good indication of the ratio between the first crown and the rosette. Only that in this case it will be indicated as 1+ratio since α2is the crown width plus α1. With this is proven that the optimal proportion is always between 0.65 and 0.75 , usually between 0.68 and 0.70 . Their similarity points toavalidity of the results. It is also observed that the dependency on α3 is little and negative, meaning that the higher α3 , the smaller α2 . The goodness of fit (whose r2 is calculated a verification of the goodness of the adjustment) proves both linear and parabolic adjustments appropriate, since it is usually close to 0.99999 and higher. Moreover, for linear adjustments, the 13th inferior value is over 0.999, and for parabolic adjustments, the 16th is over 0.9999. A few cases are shown here, the rest are presented in the appendix. Here,bcases with one crown and two crowns with4elements are presented. 3.10.1 Adjustment results for one-crown cases Since the optimal proportion for a case is not the same as the one for its two-crown counterpart (along its maximal variation and criterion 2, which we will define later) we can affirm that the second crown changes the shape of the curves even ifits precise value α3 does not affect too much. In conclusion: its presence changes the result, but its precise value does not affect it. Table 3.1 Adjustment results for one-crown cases. 4E type 1 type 0 type 2 type 0 type3type 0 type4type 0 type5type 0 p1 p1 p10.683569908 0.693775343 0.693285007 0.754764074 0.693582651 p0 p0 p0-0.014835322 -0.006490713 -0.005586549 -0.02821319 -0.005747557 6E type 1 type 0 type 2 type 0 type3type 0 type4type 0 type5type 0 p1 p1 p10.699419006 0.699419006 0.701992941 0.747145291 0.69985529 p0 p0 p0-0.00590665 -0.00590665 -0.006206542 -0.022657676 -0.005526055 8E type 1 type 0 type 2 type 0 type3type 0 type4type 0 type5type 0 p1 p1 p10.701825963 0.701825963 0.704007907 0.774798469 0.70255688 p0 p0 p0-0.005596863 -0.005596863 -0.005826572 -0.032972224 -0.005686963 12E type 1 type 0 type 2 type 0 type3type 0 type4type 0 type5type 0 p1 p1 p10.7041073 0.7041073 0.706035086 0.756101823 0.755717653 p0 p0 p0-0.005666085 -0.005666085 -0.005878402 -0.024650382 -0.024717687 24E type 1 type 0 type 2 type 0 type3type 0 type4type 0 type5type 0 p1 p1 p10.697316181 0.705272161 0.706944907 0.760888331 0.756241747 p0 p0 p0-0.002660257 -0.00558706 -0.00576177 -0.026228621 -0.024480946 66 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 3.10.2 Adjustment results for two-crown cases Table 3.2 Adjustment results for two-crown cases. 1 1 1st crown type 4E 1 2 3 4 5 1p11.7044291 1.7044291 1.7078642 1.6705268 1.6697719 p2-0.0031716 -0.0031716 -0.0030514 -0.0033376 -0.0028602 p0-0.0093358 -0.0093358 -0.0098965 -0.0001454 -0.0009216 2p11.7044291 1.7044291 1.7078642 1.6705268 1.6697881 p2-0.0031716 -0.0031716 -0.0030514 -0.0033376 -0.0027438 p0-0.0093358 -0.0093358 -0.0098965 -0.0001454 -0.0009540 2nd 2nd 2nd 3p11.7019359 1.7019359 1.7052277 1.6684431 1.6678812 crown p2-0.0005714 -0.0005714 -0.0005031 -0.0017006 -0.0011116 type p0-0.0100269 -0.0100269 -0.0105351 -0.0004649 -0.0012308 4p11.6772839 1.6772839 1.6792787 1.6829170 1.6885479 p2-0.0025825 -0.0025825 -0.0025204 -0.0024255 -0.0019403 p0-0.0014264 -0.0014264 -0.0015817 -0.0020208 -0.0043773 5p11.7076193 1.7076193 1.6788509 1.6826931 1.6880808 p20.0002524 0.0002524 -0.0020115 -0.0020918 -0.0013794 p0-0.0121905 -0.0121905 -0.0021855 -0.0027197 -0.0052554 3.11 Scores Let us defineacriterion to find the most advantageous case to use. There are two possible answers: • The stability criterion: we are searchingaconfiguration in which the effect of the crown(s) improves the result, but doesn’t dramatically affect it, since it would result in a fast deviation from the desired result. That is to say, the sensibility of the result must be very slightly dependant on the crown size, thus allowing several different sizes (without much precision in its building) to grant the right solution. Graphically, it can be seen asavery slight negative slope at the cut-off points between the surfaces/curves and the right solution surface/line. This can be spotted at the code as ’criterio_2’. It is the only criterion with a simple physical interpretation. As we can see in the figure: right solution line θ α2 uct Figure 3.42 Case with fast divergence. right solution line θ α2 uct Figure 3.43 Case with slow divergence. 3.11 Scores 67 In the figure on the left, the solution diverges rapidly whereas in the figure on the right, the simulation result remains closer to the right solution, hence, it is better to obtain a graphic like that on the right. To model this, we must first calculate the angle of the curve (its slope) in the cut-off points, then average them out, express the result as a value between 0 and 1 (that is the reason why the criterion is expressed as (π/2−θ)/(π/2), since the best case would be having a near-zero theta). To apply the same criterion to two-crown cases we can average the slopes of each cut-off point (in the code as criterio2 ). Yet this would not be a good method. Since the definition domain of the surfaces is triangular (as we can see in the previous figures), there are more cut-off points near the base. These are the cases with small α1 , which allows much more α2 - α3 configurations to appear, since they have more space to appear. This results in an over-sampling of the region near the base, which tends to have worst results (much higher slopes since they are the first cases to be considered valid, hence not very adequate) and thus artificially worsen the result, since there are more instances of the worst results. To solve this, we can firstly average the results for each α1 and then average those results, which is more comparable to the method for one crown. • The convenience criterion: in a more practical sense, we can also achieve good results with another configuration. For example, taking into account how many different α1 values are able to achieve the right solution. Also, if the regression has a better goodness of fit (for either linear or parabolic adjustments), the right-solution configurations are easier to achieve with enough adjustment precision. Finally, it is understood that using the least possible amount of nodes is preferable. These are in the code: criterio_1, criterio_3a, criterio_3b, criterio_4. The number 1131 is the highest number of nodes used in any case. The convenience criterion is be the average of any of these criteria plus the stability criterion. Criterio 1 can be appreciated in the previous images of the typologies: it is evident that the cases with more valid α1 (the typologies on the right) values are preferable over those that are only valid for few values (the typologies on the left). Criterio 3 a and b do not affect much the final result, since the values are not different among them, but can serve as indicators of the goodness of fit of the adjustment. Since they correspond to the r2 of the adjustment and are all generally over 0.9999 . They are presented in tables in the appendix. All the criteria are expressed between 0 and 1 to have an nondimensional value. We present the results with the improved criterion 2 (which seems more logical). The tables with the normal criterion 2 are analyzed, but can be consulted in the first appendix. 68 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 3.12 Table of results for a cases (improved stability criterion) The tables have a heat-map-like color code affecting all the tables in one page as a whole. Highest-scoring cases are underlined in green, and lowest-scoring cases in red. Orange is for average cases. Table 3.3 Scores for cases a with improved stability criterion. 1st Crown 4E Stab. 1 2 3 4 5 00.662867074 0.664792432 0.664880699 0.662006348 0.663763655 10.669596687 0.669597507 0.67063037 0.655767673 0.658285083 2nd Crown 20.666746609 0.669596687 0.67063037 0.655767672 0.657956754 30.665668339 0.665668339 0.669597475 0.652278565 0.653957958 40.663314223 0.663314223 0.664188574 0.664166621 0.669400217 50.679529819 0.679529819 0.664875825 0.663602237 0.669189583 1st Type 6E Stab. 1 2 3 4 5 00.663900487 0.663900487 0.664808487 0.663437268 0.663609073 10.672372249 0.672372249 0.67336265 0.660677565 0.665432649 2nd Crown 20.672372249 0.672372249 0.67336265 0.660677565 0.725558507 30.669517536 0.669517536 0.670529067 0.65737483 0.662026827 40.669042649 0.669042649 0.669803571 0.66833265 0.673140343 50.677229499 0.736839107 0.678472888 0.668293273 0.672097614 1st Type 8E Stab. 1 2 3 4 5 00.662714516 0.662714516 0.663577814 0.66264669 0.663753458 10.669474544 0.669474544 0.670253064 0.664694988 0.669413177 2nd Crown 20.669474544 0.669474544 0.670253064 0.664694988 0.669375441 30.666464867 0.666464867 0.667269408 0.661743902 0.666460208 40.668699631 0.719157603 0.669364767 0.669093893 0.669244805 50.668784391 0.668784391 0.669562523 0.668890984 0.668951643 1st Type 12E Stab. 1 2 3 4 5 00.663731231 0.663731231 0.664507975 0.663651941 0.64748072 10.669816482 0.724331233 0.670514614 0.665615869 0.671238851 2nd Crown 20.669816482 0.724331233 0.670514614 0.665615869 0.671260098 30.667016146 0.720869028 0.66774757 0.662793084 0.668798475 40.669229598 0.716622919 0.669760242 0.66970683 0.674002723 50.672713457 0.720465668 0.673489099 0.669537201 0.673512846 1st Type 24E Stab. 1 2 3 4 5 00.670563026 0.663725995 0.66448514 0.663651941 0.646842933 10.669758419 0.723107405 0.670399171 0.665200955 0.674874837 2nd Crown 20.669758419 0.723107405 0.670399171 0.665200955 0.674874055 30.667202579 0.71989872 0.667883028 0.662478884 0.672819212 40.679086512 0.712351091 0.670190515 0.669600633 0.673898285 50.669966317 0.669966317 0.671234838 0.669743844 0.673977912 3.13 Table of results for a cases (convenience criterion, parabolic adjustment, improved criterion 2) 69 3.13 Table of results for a cases (convenience criterion, parabolic adjustment, improved criterion 2) The linear adjustment is not here presented, since choosing either parabolic or linear adjustment have little final effect: both have similar (and non-decisive) values. Table 3.4 Scores for cases a with improved convenience criterion. 1st Type 4E Par. Conv. 1 2 3 4 5 00.818455561 0.807222396 0.806359551 0.805642849 0.806523585 10.749995655 0.751764206 0.751137676 0.718760524 0.719831905 2nd Crown 20.751050962 0.751763481 0.751137676 0.718760524 0.719749828 30.74989722 0.74989722 0.749995276 0.716119867 0.717423816 40.77708693 0.77708693 0.776421344 0.774647057 0.750388102 50.69824967 0.69824967 0.749257423 0.775832251 0.750329802 1st Type 6E Par. Conv. 1 2 3 4 5 00.805672047 0.805672047 0.804573127 0.804232081 0.748719416 10.778466929 0.778466929 0.777388666 0.745112673 0.774963593 2nd Crown 20.778467323 0.778467323 0.777388666 0.745112673 0.73346678 30.776427384 0.776426981 0.77535401 0.741634436 0.772343761 40.804086217 0.804086217 0.802950174 0.799930094 0.776890977 50.779239423 0.794141825 0.778224018 0.802130678 0.776630248 1st Type 8E Par. Conv. 1 2 3 4 5 00.804051489 0.804051489 0.802499006 0.802266414 0.802541549 10.748196191 0.775789999 0.746622417 0.771242931 0.69174163 2nd Crown 20.748196191 0.775789999 0.746622417 0.771242931 0.71950997 30.745675417 0.77326915 0.744108147 0.766967897 0.688792953 40.801789972 0.813423243 0.800187936 0.79658371 0.745928874 50.747582136 0.747582136 0.746008319 0.79763819 0.774516462 1st Type 12E Par. Conv. 1 2 3 4 5 00.801653383 0.801653383 0.795216277 0.796328469 0.792284657 10.772523324 0.813665967 0.770045337 0.791293084 0.739353998 2nd Crown 20.772523324 0.813665967 0.770045337 0.791293084 0.73935931 30.769170721 0.810147711 0.766701057 0.785282412 0.734765125 40.794848491 0.805806167 0.789676182 0.785458029 0.765612339 50.772805448 0.812275368 0.742569075 0.788294882 0.773005343 1st Type 24E Par. Conv. 1 2 3 4 5 00.739788738 0.79369446 0.778853309 0.796328469 0.778421109 10.706343145 0.774976441 0.701198291 0.791189356 0.76229367 2nd Crown 20.706343145 0.774976441 0.701198291 0.791189356 0.762293474 30.700399146 0.7688691 0.695264217 0.785203862 0.744538383 40.785377242 0.792817333 0.767238057 0.785431479 0.754976168 50.733509286 0.733509286 0.728521394 0.788346543 0.762511545 70 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements 3.14 Table of results for b cases (stability criterion, improved criterion 2) Table 3.5 Scores for cases b with improved stability criterion. 1st Type 4E Stab. 1 2 3 4 5 00.662867074 0.664792432 0.664880699 0.661266908 0.665107457 10.669596687 0.669597507 0.67063037 0.655125352 0.657069575 2nd Crown 20.666746609 0.669596687 0.67063037 0.655125352 0.656748736 30.665668339 0.665668339 0.669597475 0.651686089 0.652751069 40.663049668 0.663049668 0.663939203 0.662503956 0.668870601 50.680489353 0.680489353 0.664810198 0.663602237 0.669189583 1st Type 6E Stab. 1 2 3 4 5 00.665375549 0.665375549 0.666223756 0.661735864 0.662289038 10.672906982 0.672906982 0.67380076 0.659856742 0.665176447 2nd Crown 20.672906982 0.672906982 0.67380076 0.659856742 0.582605066 30.669961078 0.669961078 0.670906241 0.656584929 0.661811744 40.668389201 0.668389201 0.669169402 0.667477718 0.672320368 50.676801156 0.676801156 0.678054659 0.667583281 0.672485301 1st Type 8E Stab. 1 2 3 4 5 00.662004178 0.662004178 0.662831016 0.661149374 0.662213948 10.669062701 0.669062701 0.669844974 0.663456302 0.669295302 2nd Crown 20.669062701 0.669062701 0.669844974 0.663456302 0.376978169 30.666102136 0.666102136 0.666905746 0.660550301 0.666358192 40.667103793 0.667103793 0.66781731 0.667919264 0.668005861 50.668656835 0.668656835 0.669470164 0.668125481 0.668525596 1st Type 12E Stab. 1 2 3 4 5 00.661888992 0.661888993 0.662711356 0.661291459 0.66143681 10.669140813 0.669140813 0.669855879 0.664483703 0.668260462 2nd Crown 20.669140813 0.669140814 0.669855879 0.664483703 0.667179674 30.66637898 0.66637898 0.667120684 0.661721335 0.665178952 40.667928569 0.667719024 0.668568954 0.668519116 0.66740647 50.672300085 0.672300085 0.669745113 0.668524538 0.668126538 1st Type 24E Stab. 1 2 3 4 5 00.655531694 0.661351013 0.662128576 0.660539344 0.660601902 10.669172468 0.669172468 0.669857606 0.664374929 0.668369872 2nd Crown 20.669172468 0.669172468 0.669857606 0.664374929 0.667288907 30.666673062 0.666673062 0.667385278 0.661711122 0.665189219 40.677961468 0.668418069 0.6691373 0.667961014 0.667964108 50.66926204 0.66926204 0.669941234 0.668120592 0.668530488 3.15 Table of results for b cases (convenience criterion, parabolic adjustment, improved criterion 2) 71 3.15 Table of results for b cases (convenience criterion, parabolic adjustment, improved criterion 2) Table 3.6 Scores for cases b with improved convenience criterion. 1st Type 4E Par. Conv. 1 2 3 4 5 00.818455561 0.807222396 0.806359551 0.799047915 0.805534059 10.749995655 0.751764206 0.751137676 0.711526582 0.690865815 2nd Crown 20.751050962 0.751763481 0.751137676 0.711526582 0.690785624 30.74989722 0.74989722 0.749995276 0.708898387 0.688459903 40.769947429 0.769947429 0.769285637 0.757653135 0.743845449 50.669827596 0.669827596 0.747914776 0.775832251 0.750777521 1st Type 6E Par. Conv. 1 2 3 4 5 00.80515901 0.80515901 0.804044809 0.794523248 0.741537457 10.777274735 0.777274735 0.776171496 0.73562364 0.714037902 2nd Crown 20.777274735 0.777274735 0.776171923 0.73562364 0.688663835 30.775212002 0.775212002 0.774122033 0.732153133 0.71142837 40.794639055 0.794639055 0.793507832 0.776727871 0.767623209 50.744502214 0.744502214 0.743489332 0.787364343 0.776727197 1st Type 8E Par. Conv. 1 2 3 4 5 00.797463886 0.797463886 0.795902248 0.785313825 0.794863809 10.741682986 0.741682986 0.740110148 0.754576071 0.686407116 2nd Crown 20.741682986 0.748535329 0.740110148 0.754576071 0.638822453 30.739174488 0.739174488 0.745343499 0.750312401 0.683462404 40.785033894 0.785033891 0.783443948 0.761144179 0.729261931 50.740918956 0.740918956 0.739353939 0.780426499 0.740444039 1st Type 12E Par. Conv. 1 2 3 4 5 00.791909008 0.791909008 0.784599127 0.716973268 0.725626055 10.707736078 0.707736078 0.733040101 0.684688282 0.723578466 2nd Crown 20.707736078 0.707736078 0.733040101 0.684688282 0.7290554 30.759948657 0.759948657 0.757481562 0.678692705 0.717945134 40.771755888 0.771703117 0.766611006 0.675234604 0.685198179 50.757671175 0.757671175 0.754379931 0.705518897 0.708145656 1st Type 24E Par. Conv. 1 2 3 4 5 00.71288085 0.714335635 0.726392689 0.666608398 0.699555499 10.71098593 0.683208152 0.70585217 0.631831845 0.705922379 2nd Crown 20.71098593 0.683208152 0.70585217 0.631831845 0.705652138 30.70505604 0.70505604 0.69992905 0.620555843 0.715663613 40.729172034 0.699008335 0.683268117 0.611434532 0.6678015 50.707029177 0.707029177 0.701893951 0.647283447 0.71016237 78 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements Figure 3.52 Increasing α1 surfaces for 6E type 5 type 2 a. Figure 3.53 adjustment plane for 6E type 5 type 2 a. Figure 3.54 Increasing α1 surfaces for 8E type 2 type 4 a. Figure 3.55 adjustment plane for 8E type 2 type 4 a. Figure 3.56 Increasing α1 curves for 6E type 3 type 5 b. Figure 3.57 adjustment plane for 6E type 3 type 5 b. 3.21 Chosen configurationgraphicresults 79 Figure 3.58 Increasing α1 surfaces for 24E type 1 type 4 b, where some invalid cases warp low α1 surfaces for low α2values. Figure 3.59 adjustment plane for 24E type 1 type 4 b. Figure 3.60 Chosen : 6E type 2 type 5 a. Figure 3.61 Chosen : 4E type 1 type 0 a. Figure 3.62 Chosen : 4E type 1 type 5 free meshing in half-hexagon area. Figure 3.63 Chosen : 4E type 1 type 5 structured meshing in half-hexagon area. Figure 3.64 Chosen : 6E type 5 type2a with free meshing. Figure 3.65 Chosen : 6E type 2 type0b with free meshing. 80 Chapter 3. First Parametric Study: effect of one and two crowns surrounding the singular elements Figure 3.66 Chosen configuration: 8E type 2 type 4 a with free meshing. Figure 3.67 Chosen configuration: 8E type 2 type4awith structured meshing. Figure 3.68 Chosen configuration: 6E type 3 type 5 b with free meshing. Figure 3.69 Chosen configuration: 6E type 3 type5bwith structured meshing. 4 Second Parametric Study: optimal α1vs hexagon side-length proportion Now , with the optimal proportion among αi values, we should find out if there is an optimal value or range of values for α1 to calculate the right solution. Even when theoretically every configuration is close enough, there may be minor differences in precision. Hence, there are two possible studies: Element-Hexagon proportion keeping the hexagon side-length constant, varying α1 values. The effect of the hexagon side-length also should be taken into account (the same study should be made for several values). Hexagon-crack proportion varying the hexagon side-length in proportion to the crack. The effect of the seeds should also be taken into account (the same study should be made for several configurations). 4.1 Element-Hexagon proportion Now, as explained before, the only parameter is α1 . For one-crown cases, α2 depends directly on α1 . And for two-crown cases, α3 is mα1 (with the effect of m being studied), and α2 chosen according to the adjustment plane. The code is basically the same as before (it is added in the appendix for closer examination), but with a minor difference: all the cuts are done in the same step (to fix the numeration problem). Besides, only b cases will be slightly modified in order to achieve an even better meshing. Using a variable amount of elements between the last crown and the half-hexagon in order to make it always be meshed with elements of roughly the same size as those in the lateral and upper sides of the hexagon can be achieved by a simple formula: using (1−α3) times the number of seeds in the oblique sides of the hexagon. Upper segments will be seeded in a similar fashion. This, however, as will be seen later, is the equivalent of changing the behaviour of the mesh for some cases, and it is easier to compute them separately, and is an indicator of whether the result is a consequence of the mesh or is error-dominated. The same cannot be applied to a cases since the half-hexagon regions with only one seed are poorly fit for a structured meshing. It is worth noting that some cases plot of uct vs α1opt will exhibit an uneven curve with severe sudden steps. The curve reaches the exact value, but not in a steady, stable, and easily-predictable manner. In those cases, two solutions have been studied. The first one is imposing a structured meshing, which sometimes results in a flattening of the curve, but often stops the curve before reaching the correct value and reduces the range of validity of the result. The second option is using an inverse method to create the crown: using α2 as a parameter, choosing α3=mα2 , and α1 thanks to the adjustment. This has been proven to augment the range of validity of our result, flatten the curve, but not completely solve the problem since sometimes it is not enough to improve the whole curve. Examples of these methods are added in the following sections. It is probably due to numerical motives involving the adjustment. To keep track of the mesh-qualities, we will calculate for each case: the α1 values that result in the right solution, the maximum and minimum value of the curve, the value of m , and relevant notes. These will be analysed later on, after graphics are presented. The range of α1 in the figures is limited not only by the hexagon region, but also by the acceptable values in the adjustment. 81 82 Chapter 4. Second Parametric Study: optimal α1vs hexagon side-length proportion 4.1.1 6E type 2 type 5 a Figure 4.1 6E type 2 type 5 a (Second Study). As seen in the graphics, the parabola becomes narrower as m increases. Its maximum value descends and its minimum value ascends (it becomes slightly flatter) but not in a relevant manner. The narrowness increase is imposed by the second cut-off point, which moves more than the first. It is also relevant that small m values highly affect the result (in the first image, the correct value is not met), while for higher values it has less impact. 4.1.2 4E type 1 type 0 a Figure 4.2 4E type 1 type 0 a free and structured (Second Study) respectively. 4.1 Element-Hexagon proportion 83 For this case, there is an irregular zone which does not avoid the curve to meet the desired value, but reduces the range of α1. A structured mesh does not solve the problem, since the desired values is not met. Figure 4.3 4E type 1 type 0 a, alternative (Second Study). In this case, unlike previous tests for other configurations, the problem is not completely solved: the validity region extends considerably further than before (previously from 0.25 on, now from 0.18 on), but is not fixed. It does not change the flatness of the curve, which is an advantage over the previous two-crown case, since a flatter curve means that it always remains close to the exact solution. The max value is the one in the regularised section, not counting the irregular section prior. Yet a few validity sections seem to appear roughly between 0.15 and 0.2, where it again matches the desired solution. 4.1.3 4E type 1 type 5 b We will find analogous results. For the free case, values bellow m=2.5 were invalid. This perhaps has to do with the fact that there were few valid α1 values for this configuration. Hence, it reaches the value with a good enough flatness but only for certain conditions.Structured meshing fixes the oscillations but makes the result quickly diverge when α1 grows, as seen in the images. The flatness of the curve is not an optical illusion due to the axis values, as seen in the close-up. And the alternative method, again, makes increase the validity region but does not regularise the whole curve. free Figure 4.4 4E type 1 type 5 b free (Second Study). 84 Chapter 4. Second Parametric Study: optimal α1vs hexagon side-length proportion structured Figure 4.5 4E type 1 type 5 b structured (Second Study). Figure 4.6 Close-up of the curve.. alternative Figure 4.7 4E type 1 type 5 b alternative (Second Study). 4.1 Element-Hexagon proportion 85 4.1.4 6E type 5 type 2 a Figure 4.8 6E type 5 type 2 (Second Study). The conclusions are similar to the first case presented. Here, a structured version of the mesh has also been studied, with virtually no differences. If anything, it appears to sometimes reduce the oscillations, but they are already not very pronounced. 86 Chapter 4. Second Parametric Study: optimal α1vs hexagon side-length proportion 4.1.5 8E type 2 type 4 a With this last analysis, we can observe that a cases have been less flat than b cases, but also more smooth and stable versus m . Besides, both options: free meshing, with few in-hexagon elements, and structured, with substantially more, have similar results. Hence this case, whose validity was doubtful (since the transition was not completely appropriate), can be considered valid. Only two cases for each option are here shown for comparison. Figure 4.9 8E type 2 type 4 free (Second Study). Figure 4.10 8E type 2 type 4 structured (Second Study). 4.1.6 6E type 3 type 5 b These cases have another beneficial quality: they are the flattest curves we have found, and even intersect the right result at several points. However, the attempts to fix the oscillations result in a major improvement of the result, being one of the closest to the right solution, specially the structured version. Using the mesh options from the previous section does not improve the results either, so the problem does not have to do with the variable number of elements for the half hexagon based on the space between the last crown and the hexagon. In the appendix some graphics are presented using that method, and no improvement is found. alternative method In comparison, the curve has lost its best trait (the proximity to the exact solution) but has won manageability, since the exact cut-off points are clear. The last image is the result of applying the alternative method to a structured configuration. There does not seem to be a big difference or improvement: both structured, and alternative plus structured have similar results. As exposed here and later, its effect at best is to flatten the curve while causing convergence issues, and at worst is to lump the curve but clarify the cut-off points. 4.1 Element-Hexagon proportion 87 Figure 4.11 6E type 3 type 5 free (Second Study). Figure 4.12 6E type 3 type 5 structured (Second Study). 94 Chapter 5. Third Parametric Study: optimal Singular Element-Crack proportion Free (a) Use of #1 α1for the free case. (b) Use of #2 α1for the free case. Structured (a) Use of #1 α1for the structured case. (b) Use of #2 α1for the structured case. 5.1.6 6E type 3 type 5 b Free (a) Use of #1 α1for the free case. (b) Use of #2 α1for the structured case. 5.1 Description 95 Structured (a) Use of #1 α1for the free case. (b) Use of #2 α1for the structured case. Alternative Figure 5.10 6E type 3 type 5 b alternative process. In this case it is more relevant how structured meshing seems to flatten the curve, an effect not that evident until now. However, this comes at a cost: for small α1 values, the limitations of out configuration prevent the element to grow up to its limit (as seen in the third image). Although this improves the results, the added difficulty does not seem to make up for the result. As a test, the results are here presented for m=2 for two crowns and m=1.5for one crown. This does not change the habitual range of validity. 5.1.7 6E type 2 type 0 b Figure 5.12 Structured. In this case, the cut-off point is close to the start of the study. It is worth noting that the structured case was invalid in the beginning. no α1 could meet the desired value with the hexagon side-length of 0.1. However, now it is possible to achieve by slightly changing that value. There are a few inconveniences in structured meshes. A slight miscalculation in the proportions may make the curve loose its parabolic shape, which makes only one cut-off point to appear, as it will be shown in the next page. Results are summed in the table below, with the 96 Chapter 5. Third Parametric Study: optimal Singular Element-Crack proportion (a) Free. (b) Alternative. percentage indicating percentage of the crack length. 5.2 Table of results 97 Table 5.1 Third study summary table. α1 Best Singular Element Size (%) Max α1within 1% accuracy (%) 6E type 2 type 5 free 0.22369874 4.48928736 14.4165444 0.34219 6.84384561 16.2561366 4E type 1 type 0 free 0.25303 5.3787546 15.1773706 0.42981 8.59627852 18.1121094 4E type1type0alt 0.24103 8.59740526 18.1131118 0.42942 8.58832105 18.1058862 4E type 1 type 5 free 0.33979 4.71241692 15.8042594 4E type1type5str 0.24988 5.02865919 16.0126245 4E type1type5alt 0.18464478 3.82144632 15.28741 6E type 5 type 2 free 0.18381 3.71832262 13.6458931 8E type 2 type 4 free 0.18474 3.64200668 13.6136358 0.33255 6.65096916 15.882028 0.18290 3.42018662 13.8166996 0.33478 6.6273513 15.9984622 6E type 3 type 5 free 0.15498 3.27773902 14.0275329 0.21069 4.30452638 14.3911798 6E type3type5str 0.11197 3.57832243 12.6156338 0.24192 4.9382774 14.7730228 6E type3type5alt 0.18769 8.65937038 18.05977937 6E type 2 type 0 free 0.23957 4.79091461 14.9687985 6E type2type0str 0.24093 4.77642037 15.1484462 6E type2type0alt 0.24093 4.79091461 15.0964382 5.2 Table of results Figure 5.13 Poorly adjusted structured case, in which the parabolic shape is lost alongside its benefits. In conclusion, the results are within a 1 % accuracy range when the singular element is less than 16 % of the crack length. The results are practically accurate when the element is around 5 % of the crack length. This, althoughagood result, is slightly worse than results put forward by [ 19 ], in which a quarter point of 20 % of the crack length can still be acceptable. This, however, might be achieved by using a natural isoparametric element and not a collapsed one, as the previous results were found employing these. It is important to note that the highest-scoring configuration (4E type1type 0) also has the widest range of validity, and coincidentally, remains close enough to the solution for the most part for small α1 . The narrower interval is due to structured meshes for 6E type3type 5. The flatter curves are those of 4E type 3 type5(which even has two valid cut-off points), although for the structured mesh it hasarange limitation. They are followed by 6E type 2 type 5. In general, it will be advisable to keep the element under 10 % of the crack length and use the aforementioned configurations. 6 Conclusions, summary, further research 6.1 Conclusions As we have seen, we can conclude that the element can properly capture the displacement tendencies near the crack tip, but to ensure a correct solution, the following recommendations should be applied: • The element should preferably measure between 4% and 8% of the crack length, but a 15% of the crack length still gives solution with 1% accuracy. •One or two crowns of non-singular elements should surround the singular elements. •The first crown size should be between 0.65 and 0.7times the singular element size. •The second crown size should be between 0.88 and 1times the first crown size. • The transitioning half-hexagon can determine the solution and its length is suggested to be 10% of the crack length. •A crown configuration of the types summarised below can improve the performance of the element. The following table presents the suggested: configurations, αi , hexagon side-length, element size, adjustment coefficients and validity intervals. 6.1.1 Developed materials Handbook rules aside, during the research process, the following tools have been created: • A MATLAB script to apply boundary conditions of any kind to a square plate, a MATLAB script to read Abaqus .inp files and carry a parametric study by comparing the result of several of them. • a MATLAB script to visually present the results of such study and find an adjustment for the right element and crown(s) proportion. • A web-page to visualize meshes and visually pick nodes from it (to make sure your boundary. conditions are being applied to the desired nodes) • More than 350 Python macros to parametrically build meshes (and the MATLAB scripts to easily change the defining commands). •More than 214023 meshes studied • And several MATLAB scripts to refine the previous conclusions (finding the right hexagon-SE proportion or the right SE-crack proportion) 99 100 Chapter 6. Conclusions, summary, further research 6.2 Summary 6.2.1 Two crown cases Firstly, α1 is chosen. Then, α3 is chosen as α3=mα1 (for instance, m=2 ). Then, α2 is calculated by the following adjustment: α2=p00 +p10α1+p01α3+p20α2 1+p11α1α3+p02α2 3 And in order to simplify the resulting expression: α2=p00 + ( p10 +mp01)α1+ ( p20 +mp11 +m2p02)α2 1 α2=p0+p1α1+p2α2 1 6.2.2 One crown cases In these cases, α2is fixed and α1is calculated as follows. Evidently: p00 =p0,p10 =p1,p20 =p2. α1=p00 +p10α2+p20α2 2 6.2.3 Inverse method As seen before, choosing α2 beforehand reduces oscillations in the result. Hence, α3=m2α2 (for instance m2=1.75 ). And finally, α1=p00 +p10α2+p01α3+p20α2 2+p11α2α3+p02α2 3 , which can be similarly expressed as: α1=p00 + ( p10 +m2p01)α2+ ( p20 +m2p11 +m2 2p02)α2 2or α1=p0+p1α2+p2α2 2 The same thing applies to one-crown cases just deleting the square term. 6.2.4 Table of results The following parameters are presented in the table on the following page: • SE size: recommended size of the singular element near the crack tip (measured from the crack vertex to its base). • hex. side: hexagon side-length used in the calculations. •α1,α2,α3: their optimal value. •pi j: adjustment coefficients. •α2−α1 : size of the first crown divided by the hexagon side-length (like α coefficients are expressed). •α2−α1 α1 : size of the first crown divided by the singular-element size. Here we can see that the transition element should measure roughly 0.68 times the SE. •α3−α2 α2−α1: size of the second crown divided by the first crown size. •p1,p2,p3: coefficients of the adjustment dependant only on one parameter. •1 % acc. interval: upper limit of the 1 % accuracy interval. An Excel version has been created. It is interactive and allows the user to change m and m2 values. This table is presented for appropriate values m=2.25 and m2=1.5. The same nomenclature for coefficients is employed in each section. However, as explained here, they multiply different variables. Also, the symbol Idem is here employed to remark that the value is equal to the one in the upper box. Furthermore, an Excel file is attached to consult the values with all their decimals. A square is added to modify the m value, hence allowing to change the proportion between α values (depending on the case). Additionally, the mutual validity of normal and alternative processes is confirmed, since their results are similar enough as seen in the table. Table 6.1 Summary. 6.2 Summary 101 TWO-CROWN CASES NORMAL PROCESS SE size hex. side α1p00 p10 p01 p20 p11 p02 α2α3α2−α1 α2−α1 α1α3−α2 α3−α2 α2−α1p0 1 pp2 α1p00 p10 p01 p20 p11 p02 α2α3α2−α1 α2−α1 α1α3−α2 α3−α2 α2−α1p0 1 pp2 α1p00 p10 p01 p20 p11 p02 α2α3α2−α1 α2−α1 α1α3−α2 α3−α2 α2−α1p0 1 pp2 α1p00 p10 p01 p20 p11 p02 α2α3α2−α1 α2−α1 α1α3−α2 α3−α2 α2−α1p0 1 pp2 1% acc. interval 6E type2type 5 a free 4.4893 0.1 0.2237 0.0920 1.4444 -0.2205 0.0554 0.3029 0.1207 0.3716 0.5033 0.1479 0.6612 0.1317 0.8905 0.0920 0.9483 1.3478 14.4165 6.8438 0.1 0.3422 Idem Idem Idem Idem Idem Idem 0.5744 0.7699 0.2322 0.6785 0.1956 0.8424 Idem Idem Idem 16.2561 4E type1type 5 b free 4.7124 0.1 0.3398 -0.0017 1.6294 0.0038 0.0871 0.0329 -0.0091 0.5681 0.7645 0.2284 0.6720 0.1964 0.8600 -0.0017 1.6379 0.1152 15.8043 4E type1type 5 b str 5.0287 0.1 0.2499 Idem Idem Idem Idem Idem Idem 0.4148 0.5622 0.1649 0.6599 0.1475 0.8943 Idem Idem Idem 16.0126 6E type5type 2 a free 3.7183 0.1 0.1838 0.0537 1.4840 -0.1386 0.0837 0.2403 0.0809 0.3041 0.4136 0.1203 0.6543 0.1095 0.9104 0.0537 1.1721 1.0340 13.6459 8E type2type 4 a free 3.6420 0.1 0.1847 0.0600 1.4878 -0.1609 0.0803 0.2354 0.1012 0.3063 0.4157 0.1215 0.6580 0.1094 0.8998 0.0600 1.1259 1.1223 13.6136 6.6510 0.1 0.3326 Idem Idem Idem Idem Idem Idem 0.5585 0.7482 0.2260 0.6795 0.1897 0.8395 Idem Idem Idem 15.8820 8E type2type 4 a str 3.4202 0.1 0.1829 Idem Idem Idem Idem Idem Idem 0.3035 0.4115 0.1206 0.6592 0.1081 0.8964 Idem Idem Idem 13.8167 6.6274 0.1 0.3348 Idem Idem Idem Idem Idem Idem 0.5627 0.7533 0.2279 0.6808 0.1906 0.8360 Idem Idem Idem 15.9985 6E type3type 5 a free 3.2777 0.1 0.1550 0.0014 1.6216 0.0004 0.1420 0.0173 -0.0042 0.2567 0.3487 0.1017 0.6565 0.0920 0.9040 0.0014 1.6226 0.1596 14.0275 4.3045 0.1 0.2107 Idem Idem Idem Idem Idem Idem 0.3504 0.4741 0.1397 0.6630 0.1237 0.8855 Idem Idem Idem 14.3912 6E type3type 5 a str 3.5783 0.1 0.1120 Idem Idem Idem Idem Idem Idem 0.1851 0.2519 0.0731 0.6532 0.0668 0.9138 Idem Idem Idem 12.6156 4.9383 0.1 0.2419 Idem Idem Idem Idem Idem Idem 0.4033 0.5443 0.1614 0.6671 0.1410 0.8739 Idem Idem Idem 14.7730 ALTERNATIVE PROCESS SE size hex. side 1% acc. interval 4E type1type 5 b alt 3.8214 0.1 0.1846 0.0013 0.6120 -0.0024 -0.0174 -0.0113 0.0053 0.3066 0.4599 0.1220 0.6605 0.1533 1.2570 0.0013 0.6084 -0.0344 15.2874 6E type3type 5 b alt 8.6594 0.1 0.1877 -0.0569 0.6245 0.1714 -0.0071 -0.0406 -0.1201 0.3157 0.4735 0.1280 0.6819 0.1578 1.2332 -0.0569 0.8815 -0.3383 18.0598 ONE-CROWN CASES NORMAL PROCESS SE size hex. side 1% acc. interval 4E type1type 0 ab free 5.3788 0.1 0.2530 0.0009 0.6179 0.0000 0.1480 0.0000 0.0000 0.4197 n.a. 0.1667 0.6588 n.a. n.a. 0.0009 0.6179 0.1480 15.1774 8.5963 0.1 0.4298 Idem Idem Idem Idem Idem Idem 0.7236 n.a. 0.2938 0.6836 n.a. n.a. Idem Idem Idem 18.1121 6E type2type 0 b free 4.7909 0.1 0.2396 0.0009 0.6256 0.0000 0.1475 0.0000 0.0000 0.3988 n.a. 0.1592 0.6646 n.a. n.a. 0.0009 0.6256 0.1475 14.9688 6E type2type 0 b str 4.7764 0.1 0.2409 Idem Idem Idem Idem Idem Idem 0.4011 n.a. 0.1602 0.6647 n.a. n.a. Idem Idem Idem 15.1484 ALTERNATIVE PROCESS SE size hex. side 1% acc. interval 4E type1type 0 ab alt 8.5974 0.1 0.2410 -0.0003 0.6162 0.0000 -0.0301 0.0000 0.0000 0.3995 n.a. 0.1584 0.6573 n.a. n.a. -0.0003 0.6162 -0.0301 18.1131 8.5883 0.1 0.4294 Idem Idem Idem Idem Idem Idem 0.7229 n.a. 0.2935 0.6835 n.a. n.a. -0.0003 0.6162 -0.0301 18.1059 6E type2type 0 b alt 4.7909 0.1 0.2409 -0.0004 0.6137 0.0000 -0.0301 0.0000 0.0000 0.4011 n.a. 0.1602 0.6648 n.a. n.a. -0.0004 0.6137 -0.0301 15.0964 102 Chapter 6. Conclusions, summary, further research 6.3 Further research The next steps regarding the research are suggested to be the following: • Convergence keeping the suggestions while augmenting the number of elements in the plate. For this, special care to the relation between the meshing parts and the size of the element should be given. With this, an improvement for small SE values could be achieved. • Study on the effect of transition region. As seen before, structured meshing affects the optimal values but can make the behaviour versus SE size more manageable. Finding the best way to transition (fixating the number of elements on the half-hexagon and its arrangement) may make the study results more conclusive. • Study on the effect of using singular-elements as transition elements in the crown. As referenced before, the dominance zone is neither singular or non-singular, hence using linear or quadratic elements may help capture the displacements, specially for high-precision studies with fine meshes. •Changing the element to give a quadratic behaviour on the opposite side instead of a linear one. • Implementing an Abaqus sub-routine to use the element directly on Abaqus, in order to implement this whole study method in one program, benefiting from its time-cost-improved algorithm. • Using an isoparametric element instead of a collapsed one, revisiting the element equations but trying to build it in another way. 6.4 Personal conclusions Thanks to this project, I feel I have acquired several new abilities. For instance, planning has proven to be a major issue, having to employ several computers at a time to reduce total time computing cases, or leaving the computer calculating during the night instead of working simultaneously on two things while keeping an eye on simulations. Furthermore, while managing a gargantuan amount of information, I learnt how important it is to keep a code for identification and sticking to it. Similarly, I learnt that using a lot of information also require useful tools to analyse it. Working on and developing such tools is fundamental and should not be avoided by trying to re-adapt other tools or doing the work manually. Finally, I learnt the importance of revising the results and keeping up to date the final product (in this case, this document). Losing the main purpose of the study (putting forward conclusions) can lead to self-indulgence in minor details. That is, perfectionism is fine as long as it does not make you lose time and forget your goals. All in all, it is not only the amount of work put into a research, but also the way everything is explained and the need of those results. 7 Appendix: Auxiliary code During the preparation of the previous studies, a graphical tool to visually select nodes in the meshes was developed. This was programmed in order to make sure that the lists of nodes generated by our code, which automatically imposes boundary conditions, corresponds to the desired selection of nodes. Although it is not necessary for the results, it is a tool that might be of help for future modifications of the code to adapt it to other cases. In the following sections, the tool’s code and purpose is explained. Hyperlinks redirect to web pages in the digital edition. It was originally developed in Spanish, so text in the images (in Spanish) will be referred to in English (as the previous sections). 7.1 Introduction In order to study the meshes, a little web application has been developped. This application is located in a server, a web hosting with the domain: www.ayrsoftsoluciones.com. The app is run from a navigator window with the following url: http://www.ayrsoftsoluciones.com/nodos. This will redirect to the screen in the left. To access, the username and password previously assigned by the administrator must be entered. It will display the right screen capture. There is a main menu with three options: 1. Administration. From this option we can register different users which will be able to use the app. 2. Structural Study. From this option we will work the different meshes. 3. Exit. To abandon the web application. (a) Access screen. (b) First screen. 7.2 Application functioning First of all, we must create a study. This study will have a name assigned and will have two archives associated containing information related to each mesh. When we click on Structural Study, we will find a screen like the left one. 103 110 Chapter 8. Appendix: Tables 1st Crown 12E Stab. 1 2 3 4 5 00.663731231 0.663731231 0.664507975 0.663651941 0.64748072 10.661010367 0.663773113 0.6615218 0.655266717 0.657883186 2nd Crown 20.661010367 0.663773113 0.6615218 0.655266717 0.657873389 30.65727861 0.660009566 0.657824262 0.651583365 0.654220569 40.657273219 0.729622618 0.65750124 0.656943804 0.656345757 50.661164036 0.663845391 0.661450194 0.656866707 0.655981296 1st Crown 24E Stab. 1 2 3 4 5 00.670563026 0.663725995 0.66448514 0.663651941 0.646842933 10.660714185 0.66341278 0.661182001 0.655266717 0.65809449 2nd Crown 20.660714185 0.66341278 0.661182001 0.655266717 0.658093391 30.657220132 0.659887725 0.657724405 0.651583365 0.654458495 40.670190027 0.730816512 0.65700697 0.656943804 0.656383712 50.660992139 0.660992139 0.661663657 0.656866707 0.656177615 8.2 Score Table: a cases Convenience (normal criterion 2) Table 8.2 Scores for a cases Convenience criterion (normal criterion 2). 1st Crown 4E Par. Conv. 1 2 3 4 5 00.818455561 0.807222396 0.806359551 0.805642849 0.806523585 10.747594263 0.749362988 0.748679328 0.716603627 0.718064673 2nd Crown 20.749362459 0.749362089 0.748679328 0.716603627 0.717982834 30.747226678 0.747226678 0.746551762 0.713772378 0.715482275 40.773869079 0.773869079 0.77311785 0.772106039 0.748031671 50.69688269 0.69688269 0.746705577 0.77331063 0.741125443 1st Crown 6E Par. Conv. 1 2 3 4 5 00.805672047 0.805672047 0.804573127 0.804232081 0.748719416 10.775594144 0.775594144 0.774418676 0.742889835 0.772079789 2nd Crown 20.775594537 0.775594537 0.774418676 0.742889835 0.742210059 30.773261304 0.773260901 0.773610233 0.739220478 0.769203924 40.800928955 0.800928955 0.799688303 0.796918959 0.773944384 50.77675642 0.77675642 0.775624608 0.799179842 0.77096719 8.3 Score Table: b cases Stability (normal criterion 2) 111 1st Crown 8E Par. Conv. 1 2 3 4 5 00.804051489 0.804051489 0.802499006 0.802266414 0.802541549 10.745791341 0.773766583 0.744163956 0.768473117 0.687408497 2nd Crown 20.745791341 0.773766583 0.744163956 0.768473117 0.715179183 30.742979209 0.770949942 0.741655322 0.763957496 0.684544826 40.798883357 0.819977597 0.797187852 0.79345364 0.743750804 50.74717619 0.74717619 0.745578814 0.794504445 0.77146372 1st Crown 12E Par. Conv. 1 2 3 4 5 00.801653383 0.801653383 0.795216277 0.796328469 0.792284657 10.770321795 0.798526437 0.767797133 0.788705796 0.736015082 2nd Crown 20.770321795 0.798526437 0.767797133 0.788705796 0.736012633 30.766736337 0.794932846 0.76422023 0.782479982 0.731120648 40.791859397 0.809056092 0.786611432 0.782267272 0.761198098 50.769918093 0.798120299 0.739559349 0.785127259 0.768622456 1st Crown 24E Par. Conv. 1 2 3 4 5 00.739788738 0.79369446 0.778853309 0.796328469 0.778421109 10.704082087 0.760052784 0.698893999 0.788705796 0.758098583 2nd Crown 20.704082087 0.760052784 0.698893999 0.788705796 0.758098308 30.697903534 0.753866351 0.692724561 0.782479982 0.739948204 40.78315312 0.797433689 0.763942171 0.782267272 0.750597525 50.731265742 0.731265742 0.726128598 0.785127259 0.758061471 8.3 Score Table: b cases Stability (normal criterion 2) Table 8.3 Scores for b cases Convenience criterion (normal criterion 2). 1st Crown 4E Stab. 1 2 3 4 5 00.662867074 0.664792432 0.664880699 0.661266908 0.665107457 10.674264335 0.674264335 0.675507538 0.64624417 0.639912003 2nd Crown 20.674264335 0.674264335 0.675507538 0.646244169 0.639716321 30.672341151 0.672341151 0.673588843 0.642046482 0.636444725 40.650125896 0.650125896 0.650679832 0.652028704 0.659195358 50.674559229 0.674559229 0.654362506 0.653515753 0.659813078 112 Chapter 8. Appendix: Tables 1st Crown 6E Stab. 1 2 3 4 5 00.665375549 0.665375549 0.666223756 0.661735864 0.662289038 10.666021054 0.666021054 0.6668002 0.650739862 0.643650612 2nd Crown 20.666021054 0.666021054 0.6668002 0.650739862 0.580380612 30.662405569 0.662405569 0.663211672 0.646710425 0.640608352 40.65575344 0.65575344 0.656218504 0.655297786 0.66031053 50.666922595 0.666922595 0.667750216 0.655591078 0.660759765 1st Crown 8E Stab. 1 2 3 4 5 00.662004178 0.662004178 0.662831016 0.661149374 0.662213948 10.659257601 0.659257601 0.659857471 0.652002474 0.65167463 2nd Crown 20.659257601 0.659257601 0.659857471 0.652002474 0.577445738 30.655152473 0.655152473 0.655770069 0.648143501 0.64909434 40.65500209 0.655266321 0.655410151 0.655040801 0.659215701 50.660829554 0.660829554 0.661484663 0.655330658 0.659829776 1st Crown 12E Stab. 1 2 3 4 5 00.661888992 0.661888993 0.662711356 0.661291459 0.66143681 10.660296948 0.660296948 0.660850756 0.65379232 0.655626391 2nd Crown 20.660296948 0.660296948 0.660850756 0.65379232 0.655372469 30.656616962 0.656616962 0.657193442 0.650180929 0.652058427 40.655629587 0.661692331 0.656036148 0.655600587 0.654003501 50.660528319 0.660528319 0.660874644 0.655645975 0.654024527 1st Crown 24E Stab. 1 2 3 4 5 00.655531694 0.661351013 0.662128576 0.660539344 0.660601902 10.660188144 0.660188144 0.660723744 0.65361888 0.655230592 2nd Crown 20.660188144 0.660188144 0.660723744 0.65361888 0.655229492 30.656762103 0.656762103 0.657318235 0.650091295 0.651802182 40.668706197 0.656033885 0.655803807 0.654867356 0.65340676 50.660242377 0.660242377 0.660773776 0.655106776 0.653468682 8.4 Score Table: b cases Convenience (normal criterion 2) 113 8.4 Score Table: b cases Convenience (normal criterion 2) Table 8.4 Scores for b cases Convenience Criterion (normal criterion 2). 1st Crown 4E Par. Conv. 1 2 3 4 5 00.755957352 0.744723522 0.743861292 0.699047915 0.705534059 10.629915079 0.629915079 0.629341707 0.69263962 0.686576422 2nd Crown 20.629915079 0.629915079 0.629341707 0.692639619 0.68652752 30.62855011 0.62855011 0.62797786 0.689821819 0.684383317 40.726312445 0.726312445 0.725566754 0.675669243 0.692815527 50.668345065 0.668345065 0.700858409 0.693945551 0.699822284 1st Crown 6E Par. Conv. 1 2 3 4 5 00.715873296 0.715873296 0.714759094 0.685148248 0.687718013 10.66696701 0.66696701 0.665835497 0.67351536 0.678353413 2nd Crown 20.666967011 0.666967011 0.665835497 0.67351536 0.657804691 30.66473688 0.66473688 0.663612103 0.669855447 0.675824491 40.702194401 0.702194401 0.700984393 0.673682888 0.696244682 50.693421462 0.693421462 0.69230211 0.684366292 0.705419745 1st Crown 8E Par. Conv. 1 2 3 4 5 00.697463886 0.697463886 0.695902248 0.675938825 0.685488809 10.6906206 0.6906206 0.689002161 0.672347534 0.668113059 2nd Crown 20.6906206 0.6906206 0.689002161 0.672347534 0.672272679 30.687825961 0.687825961 0.686211953 0.667845621 0.665257552 40.682008468 0.682074523 0.680342159 0.648549563 0.657619947 50.679133076 0.679133076 0.677528505 0.667852794 0.66882564 1st Crown 12E Par. Conv. 1 2 3 4 5 00.692124897 0.692124897 0.685678127 0.686804457 0.682702052 10.701929582 0.721495855 0.699404114 0.672017178 0.656945361 2nd Crown 20.701929582 0.721495855 0.699404114 0.672017178 0.656942911 30.698343752 0.71790069 0.695826816 0.665791515 0.65205058 40.682459273 0.699403505 0.677210535 0.672866311 0.643135169 50.690534919 0.708731186 0.687954559 0.675726604 0.66002896 114 Chapter 8. Appendix: Tables 1st Crown 24E Par. Conv. 1 2 3 4 5 00.669171491 0.684168966 0.669318068 0.686804457 0.668887595 10.662398985 0.680584638 0.657210186 0.672017178 0.660866625 2nd Crown 20.662398985 0.680584638 0.657210186 0.672017178 0.66086635 30.656219929 0.674396512 0.651040211 0.665791515 0.642715866 40.673765072 0.687850213 0.654540095 0.672866311 0.642004397 50.661805036 0.661805036 0.656667205 0.675726604 0.649468373 8.5 Adjustment coefficients for a cases 115 8.5 Adjustment coefficients for a cases Table 8.5 One-crown a cases. 4E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.699316396 0.721370467 0.725865491 0.693738188 0.692766707 p p p0 0 0-0.008376937 -0.016454341 -0.017262345 -0.005477516 -0.005824042 6E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.714241418 0.714241418 0.71739044 0.700172979 0.700123043 p p p0 0 0-0.011332475 -0.011332475 -0.011799581 -0.005750364 -0.006089905 8E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.709116282 0.709116282 0.711741058 0.702693161 0.702223561 p p p0 0 0-0.00849959 -0.00849959 -0.008855271 -0.005964584 -0.00598857 12E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.711496108 0.711496108 0.713803468 0.704545669 0.703920464 p p p0 0 0-0.008575322 -0.008575322 -0.008900085 -0.006041095 -0.006002726 24E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.709116282 0.709116282 0.711741058 0.702693161 0.702223561 p p p0 0 0-0.00849959 -0.00849959 -0.008855271 -0.005964584 -0.00598857 Table 8.6 Two crown a cases. 1st crown type 4E 1 2 3 4 5 1p11.6867245 1.6867244 1.6894049 1.6703374 1.6697376 p2-0.0017726 -0.0017726 -0.0017083 -0.0032327 -0.0031663 p0-0.0050516 -0.0050517 -0.0053385 -0.0005834 -0.0008996 2p11.6871925 1.6867245 1.6894049 1.6703374 1.6697527 p20.0003820 -0.0017726 -0.0017083 -0.0032327 -0.0030488 p0-0.0057496 -0.0050516 -0.0053385 -0.0005834 -0.0009332 2nd 3p11.6846197 1.6846197 1.6867245 1.6682768 1.6678813 crown p20.0003312 0.0003312 -0.0017726 -0.0015864 -0.0014140 type p0-0.0054791 -0.0054791 -0.0050516 -0.0009105 -0.0012117 4p11.6769952 1.6769952 1.6791437 1.6827517 1.6883567 p2-0.0025177 -0.0025177 -0.0024742 -0.0023952 -0.0018377 p0-0.0018227 -0.0018227 -0.0019846 -0.0024714 -0.0048130 5p11.7074728 1.7074728 1.6787648 1.6826931 1.6880808 p20.0006530 0.0006530 -0.0021046 -0.0020918 -0.0013794 p0-0.0126425 -0.0126425 -0.0023038 -0.0027197 -0.0052554 116 Chapter 8. Appendix: Tables 1st crown type 6E 1 2 3 4 5 1p11.6961270 1.6961271 1.6976988 1.6826332 1.6849768 p2-0.0010747 -0.0010750 -0.0010529 -0.0026601 -0.0024295 p0-0.0060096 -0.0060093 -0.0061239 -0.0018505 -0.0024697 2p11.6961271 1.6961271 1.6976988 1.6826332 1.6835046 p2-0.0010750 -0.0010750 -0.0010529 -0.0026601 -0.0196051 p0-0.0060093 -0.0060093 -0.0061239 -0.0018505 0.0093866 2nd 3p11.6938452 1.6938452 1.6953310 1.6803116 1.6828005 crown p20.0006449 0.0006449 0.0006574 -0.0010668 -0.0008077 type p0-0.0062827 -0.0062827 -0.0063859 -0.0021342 -0.0027296 4p11.6901831 1.6901831 1.6917266 1.6913863 1.6968996 p2-0.0022881 -0.0022880 -0.0022329 -0.0018295 -0.0014130 p0-0.0031221 -0.0031221 -0.0032433 -0.0033897 -0.0056025 5p11.7062876 1.7005141 1.7082988 1.6913545 1.6970652 p2-0.0006975 -0.0218790 -0.0006694 -0.0017554 -0.0009870 p0-0.0087947 0.0075208 -0.0090019 -0.0034818 -0.0060645 1st crown type 8E 1 2 3 4 5 1p11.6960794 1.6960793 1.6979643 1.6882044 1.6966203 p2-0.0018152 -0.0018103 -0.0018200 -0.0027666 -0.0019847 p0-0.0050894 -0.0050928 -0.0052551 -0.0022779 -0.0050105 2p11.6960794 1.6960793 1.6979643 1.6882044 1.6965937 p2-0.0018152 -0.0018103 -0.0018200 -0.0027666 -0.0019462 p0-0.0050894 -0.0050928 -0.0052551 -0.0022779 -0.0050268 2nd 3p11.6936084 1.6936083 1.6953341 1.6858175 1.6940070 crown p2-0.0001006 -0.0000957 -0.0000581 -0.0011423 -0.0002780 type p0-0.0053612 -0.0053646 -0.0055433 -0.0025743 -0.0052755 4p11.6924175 1.6885851 1.6943677 1.6932596 1.6986479 p2-0.0019520 -0.0193744 -0.0018183 -0.0017201 -0.0014134 p0-0.0033424 0.0089982 -0.0035742 -0.0034970 -0.0055818 5p11.6981914 1.6981914 1.6996977 1.6932205 1.6985796 p2-0.0011186 -0.0011186 -0.0010269 -0.0016343 -0.0013586 p0-0.0057929 -0.0057929 -0.0059644 -0.0035734 -0.0056695 8.5 Adjustment coefficients for a cases 117 1st crown type 12E 1 2 3 4 5 1p11.6999712 1.6988288 1.7015616 1.6913745 1.7188120 p2-0.0013853 -0.0124948 -0.0013231 -0.0022309 -0.0014770 p0-0.0056876 0.0021801 -0.0058687 -0.0028155 -0.0109862 2p11.6999712 1.6988288 1.7015616 1.6913745 1.7188125 p2-0.0013853 -0.0124948 -0.0013231 -0.0022309 -0.0014512 p0-0.0056876 0.0021801 -0.0058687 -0.0028155 -0.0110062 2nd 3p11.6974904 1.6963475 1.6990279 1.6891046 1.7167820 crown p20.0002027 -0.0109115 0.0002470 -0.0007732 -0.0001593 type p0-0.0059221 0.0019489 -0.0060898 -0.0030457 -0.0111703 4p11.6940532 1.6981477 1.6955171 1.6947369 1.7188618 p2-0.0017718 -0.0242463 -0.0017068 -0.0018564 -0.0015532 p0-0.0034616 0.0102150 -0.0036098 -0.0034586 -0.0108669 5p11.7015837 1.7008330 1.7020216 1.6947319 1.7186351 p2-0.0011679 -0.0125103 -0.0011583 -0.0018965 -0.0012491 p0-0.0060627 0.0017775 -0.0060170 -0.0034501 -0.0110586 1st crown type 24E 1 2 3 4 5 1p11.7009137 1.6995978 1.7023247 1.6919146 1.7235422 p2-0.0012056 -0.0123571 -0.0011483 -0.0021321 -0.0012502 p0-0.0057354 0.0020330 -0.0058955 -0.0028287 -0.0120399 2p11.7009137 1.6995978 1.7023247 1.6919146 1.7235411 p2-0.0012056 -0.0123571 -0.0011483 -0.0021321 -0.0012496 p0-0.0057354 0.0020330 -0.0058955 -0.0028287 -0.0120399 2nd 3p11.6985027 1.6971862 1.6998871 1.6896712 1.7219162 crown p20.0002696 -0.0108863 0.0003058 -0.0007131 -0.0002321 type p0-0.0059396 0.0018318 -0.0060865 -0.0030453 -0.0121463 4p11.7146987 1.6881416 1.6967477 1.6955540 1.7198921 p2-0.0011859 -0.0124858 -0.0017459 -0.0018422 -0.0015068 p0-0.0092715 0.0046742 -0.0036277 -0.0034732 -0.0109356 5p11.7010286 1.7010286 1.7025776 1.6955513 1.7198391 p2-0.0013728 -0.0013728 -0.0014492 -0.0018151 -0.0014480 p0-0.0056455 -0.0056455 -0.0057185 -0.0035038 -0.0109780 118 Chapter 8. Appendix: Tables Table 8.7 One-crown b cases. 4E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.683569908 0.693775343 0.693285007 0.754764074 0.693582651 p p p0 0 0-0.014835322 -0.006490713 -0.005586549 -0.02821319 -0.005747557 6E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.699419006 0.699419006 0.701992941 0.747145291 0.69985529 p p p0 0 0-0.00590665 -0.00590665 -0.006206542 -0.022657676 -0.005526055 8E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.701825963 0.701825963 0.704007907 0.774798469 0.70255688 p p p0 0 0-0.005596863 -0.005596863 -0.005826572 -0.032972224 -0.005686963 12E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.7041073 0.7041073 0.706035086 0.756101823 0.755717653 p p p0 0 0-0.005666085 -0.005666085 -0.005878402 -0.024650382 -0.024717687 24E type 1 type 0 type 2 type 0 type 3 type 0 type 4 type 0 type 5 type 0 p p p1 1 10.697316181 0.705272161 0.706944907 0.760888331 0.756241747 p p p0 0 0-0.002660257 -0.00558706 -0.00576177 -0.026228621 -0.024480946 8.6 Adjustment coefficients for b cases Table 8.8 Two-crown b cases. 1st crown type 4E 1 2 3 4 5 1p11.7044291 1.7044291 1.7078642 1.6705268 1.6697719 p2-0.0031716 -0.0031716 -0.0030514 -0.0033376 -0.0028602 p0-0.0093358 -0.0093358 -0.0098965 -0.0001454 -0.0009216 2p11.7044291 1.7044291 1.7078642 1.6705268 1.6697881 p2-0.0031716 -0.0031716 -0.0030514 -0.0033376 -0.0027438 p0-0.0093358 -0.0093358 -0.0098965 -0.0001454 -0.0009540 2nd 3p11.7019359 1.7019359 1.7052277 1.6684431 1.6678812 crown p2-0.0005714 -0.0005714 -0.0005031 -0.0017006 -0.0011116 type p0-0.0100269 -0.0100269 -0.0105351 -0.0004649 -0.0012308 4p11.6772839 1.6772839 1.6792787 1.6829170 1.6885479 p2-0.0025825 -0.0025825 -0.0025204 -0.0024255 -0.0019403 p0-0.0014264 -0.0014264 -0.0015817 -0.0020208 -0.0043773 5p11.7076193 1.7076193 1.6788509 1.6826931 1.6880808 p20.0002524 0.0002524 -0.0020115 -0.0020918 -0.0013794 p0-0.0121905 -0.0121905 -0.0021855 -0.0027197 -0.0052554 8.6 Adjustment coefficients for b cases 119 1st crown type 4E 1 2 3 4 5 1p11.6963644 1.6963644 1.6979081 1.6827237 1.6851085 p2-0.0014079 -0.0014079 -0.0013688 -0.0027481 -0.0027112 p0-0.0056197 -0.0056197 -0.0057414 -0.0013880 -0.0019108 2p11.6963644 1.6963644 1.6979081 1.6827237 1.4226718 p2-0.0014079 -0.0014079 -0.0013688 -0.0027481 0.1696187 p0-0.0056197 -0.0056197 -0.0057414 -0.0013880 -0.0458324 2nd 3p11.6934941 1.6934941 1.6955364 1.6803832 1.6829188 crown p20.0003033 0.0003033 0.0003573 -0.0011741 -0.0011093 type p0-0.0057921 -0.0057921 -0.0060162 -0.0016572 -0.0021583 4p11.6903826 1.6903826 1.6919180 1.6912877 1.6973923 p2-0.0023728 -0.0023728 -0.0023251 -0.0020264 -0.0012877 p0-0.0026747 -0.0026747 -0.0027890 -0.0027916 -0.0053259 5p11.7065202 1.7065202 1.7084621 1.6912537 1.6966968 p2-0.0009088 -0.0009088 -0.0008964 -0.0020129 -0.0012330 p0-0.0083248 -0.0083248 -0.0085068 -0.0028256 -0.0056375 1st crown type 8E 1 2 3 4 5 1p11.6964109 1.6964109 1.6981849 1.6884181 1.6972325 p2-0.0019570 -0.0019570 -0.0019055 -0.0028454 -0.0022469 p0-0.0046484 -0.0046484 -0.0048355 -0.0018030 -0.0045800 2p11.6964109 1.6964109 1.6981849 1.6884181 1.2072862 p2-0.0019570 -0.0019570 -0.0019055 -0.0028454 0.3958576 p0-0.0046484 -0.0046484 -0.0048355 -0.0018030 -0.1320626 2nd 3p11.6938288 1.6938288 1.6955355 1.6860232 1.6946064 crown p2-0.0002127 -0.0002127 -0.0001703 -0.0012595 -0.0005486 type p0-0.0049281 -0.0049281 -0.0051060 -0.0020750 -0.0048363 4p11.6929775 1.6932707 1.6945762 1.6934519 1.6988877 p2-0.0019315 -0.0021166 -0.0018921 -0.0019374 -0.0015389 p0-0.0029760 -0.0029289 -0.0031070 -0.0029112 -0.0050890 5p11.6983308 1.6983308 1.6998433 1.6934486 1.6987667 p2-0.0015662 -0.0015662 -0.0014814 -0.0018590 -0.0014089 p0-0.0051553 -0.0051553 -0.0053219 -0.0029589 -0.0052688 126 Chapter 9. Appendix: Codes 6nofcrowns=2; 7else 8nofcrowns=1; 9end 10 tipo1=str2num(id( separadores (2)+1:separadores (3)−1)); 11 tipo2=str2num(id( separadores (4)+1:end)); Script 9.2 Code to define cut vectors. 9.1.1 step_03_name_archives 1clear all ,close all , 2% clc 3%% Descripcion 4%resuelve la placa rectangular inicial con grieta 5id_real_index = strvcat ( 6E_tipo_2_tipo_5 ); 6nota_index =[]; 7t_index =[]; 8criterio_1_index =[]; 9criterio_2_index =[]; 10 criterio_3a_index =[]; 11 criterio_3b_index =[]; 12 criterio_4_index =[]; 13 14 %% 15 16 for puntero=1:size ( id_real_index ,1) 17 clearvars −except t_index puntero id_real_index nota_index criterio_1_index criterio_2_index criterio_3a_index criterio_3b_index criterio_4_index ; close all ; 18 %% Identificaci ón de la carpeta que analizar 19 id=id_real_index (puntero ,:) ; %% Nombre de la carpeta donde están los .inp 20 id 21 22 %% Datos a introducir 23 hex_sidelength =0.1; % Introducir lado del hexágono 24 typeelement=1; % 1: con elementos singulares // 0: sin elementos singulares 25 vct=[0.00 0.00 0.00]; % vector con las coordenadas del vé rtice de la grieta 26 G=1; k_robin=1; NUM=400; 27 rcts =2.25596; % Right crack tip solution 28 29 %% Inicializar variables 30 [name_index, file_index ,alpha1_index,alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , singular_FE_sidelength_index ,new_entrance_inp_index] = step_01_set_variables () ; 31 32 %% Obtener parámetros del nombre 33 [nofelements , nofcrowns, tipo1 , tipo2 , separadores ] = step_02_obtain_parameters (id) ; 34 35 % cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros \ , id (1: separadores (1)−1) , \ , id ]) ; 36 % load([id , . mat ]) ; 37 38 %% Índice de nombres 39 [alpha1_index,alpha2_index,alpha3_index, singular_FE_sidelength_index , file_index ,name_index,ncasos] = step_03_name_archives(nofelements ,nofcrowns,tipo1 , tipo2 , hex_sidelength , id ,name_index, file_index , singular_FE_sidelength_index ,alpha1_index,alpha2_index,alpha3_index, separadores ) ; 40 41 % Bucle para aná lisis de todos los casos de un estudio 42 for s=1:size ( file_index ,1) 43 close all ; 44 45 %% Cambio de directorio de trabajo (cambiar al del ordenador que lo utilice ) 46 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 47 48 %% Particularizamos los valores de ciertos parámetros al caso del índice estudiado 49 alpha1=alpha1_index (:, s); 50 SE_sidelength=singular_FE_sidelength_index (s); 51 file = file_index (s ,:) ; 52 inp_malla= fileread ( file ); 53 tic 9.1 Codes for study 1 127 54 55 %% Obtenemos mallas 56 cd([ 'C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\Codigo Matlab']) ; 57 [coordenadas,coordenadasN5,coordenadas5, triangulos , triangulos5 ]=step_04_mesher(inp_malla,nofelements,SE_sidelength, alpha1,typeelement, vct); 58 59 %% Imponer condiciones de contorno 60 [neumann,neumann5,dirichlet, dirichlet5 , robin ,robin5,nodos_calculo]=step_05_boundary_conditions(coordenadas,vct , triangulos , triangulos5 ); 61 62 %% Contar elementos 63 [numNodos,numTriangulos,numTriangulos5,numLadosNeu,numLadosNeu5,numLadosDir,numLadosDir5,numLadosRob, numLadosRob5]=step_06_count_elements(coordenadas,triangulos,triangulos5,neumann,neumann5,dirichlet,dirichlet5, robin,robin5); 64 65 %% Matriz K_mn del elemento singular 66 RobMat=step_07_robin_matrix(); 67 68 %% Inicializacion de variables 69 [K,b,KROBIN,u] = step_08_initialize_system_matrixes (numNodos); 70 71 %% Generacion de la matriz de rigidez del sistema 72 K = step_09_stiffness_matrix (coordenadas, triangulos , triangulos5 ,numTriangulos,numTriangulos5,G,K); 73 74 %% Generacion de las matrices que hay que ensamblar con la matriz de rigidez y que son debida a las condiciones robin 75 KROBIN = step_10_robin_stiffness_matrix(numLadosRob,numLadosRob5,robin,robin5,coordenadas,k_robin,RobMat, numTriangulos,triangulos,KROBIN,NUM); 76 77 %% ensamblaje de las condiciones robin en la matriz de rigidez : 78 K=K+KROBIN; 79 80 %% Condiciones neumann 81 b = step_11_neumann_conditions(b,numLadosNeu,numLadosNeu5,neumann,neumann5,coordenadas,numTriangulos,RobMat, numTriangulos5,triangulos,triangulos5); 82 83 %% condiciones dirichlet 84 [b,Dir,Dir5] = step_12_dirichlet_conditions (b,numLadosDir,numLadosDir5,dirichlet, dirichlet5 ,coordenadas, numTriangulos,RobMat,numTriangulos5,triangulos, triangulos5 ); 85 86 %% vector independiente 87 b=b−K∗u; 88 89 %% solucion 90 nodosDirTotales=[Dir;Dir5]; %forma una matriz con todos los nodos con condiciones Dirichlet 91 nodosLibres= setdiff (1:numNodos,nodosDirTotales); %excluye a los nodos con condiciones Dirichlet 92 u(nodosLibres)=K(nodosLibres,nodosLibres)\b(nodosLibres); 93 t=toc; 94 t_index=[t_index ; t ]; 95 96 %% desplazamientos entorno a la grieta 97 desplazamientos_calculo =[ desplazamientos_calculo ; u(nodos_calculo) ']; 98 coordenadas_calculo=[coordenadas_calculo ; coordenadas(nodos_calculo ,1) ']; 99 end 100 101 %% Ajuste 102 cd([ 'C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\Codigo Matlab']) ; 103 [ ajuste , ajuste_2 , slope_index,max_2nd_crown_change,casos_buenos,nota_a,nota_b, criterio_1 , criterio_2 , criterio_3a , criterio_3b , criterio_4 ,criterio_1_num , criterio_1_den , criterio_5 , Coeff_vectora , Coeff_vectorb ,mean_slopes_index , criterio_2_mejorado ] = step_13_numerical_regressions (alpha1_index,alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , tipo1 , tipo2 ,id , separadores ,coordenadas,2.25596, t) ; 104 105 Coeff_vectora 106 Coeff_vectorb 107 108 109 110 if nofcrowns==2, 111 strvcat ([ 'p00='num2str(Coeff_vectorb .p00,50) ],[ 'p10='num2str(Coeff_vectorb .p10,50) ],[ 'p01='num2str( Coeff_vectorb .p01,50) ],[ 'p20='num2str(Coeff_vectorb .p20,50) ],[ 'p11='num2str(Coeff_vectorb .p11,50) ],[ 'p02=' num2str(Coeff_vectorb .p02,50)]) 128 Chapter 9. Appendix: Codes 112 else 113 strvcat ([ p1= num2str(ajuste_2 .p1,50) ],[ p2= num2str(ajuste_2 .p2,50) ],[ p3= num2str(ajuste_2 .p3,50) ]) 114 strvcat ([ p1= num2str( ajuste .p1,50) ],[ p2= num2str( ajuste .p2,50) ],[ p3= num2str( ajuste .p3,50) ]) 115 end 116 117 118 nota_index=[nota_index; nota_a, nota_b]; 119 120 %% Guardado 121 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 122 save(id); 123 cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\Codigo Matlab ]) ; 124 end Script 9.3 Code to name and identify archives from Abaqus. 9.1.2 step 01 1function [name_index, file_index ,alpha1_index,alpha2_index,alpha3_index, desplazamientos_calculo , coordenadas_calculo , singular_FE_sidelength_index ,new_entrance_inp_index] = step_01_set_variables () 2name_index=[]; 3file_index =[]; 4alpha1_index =[]; 5alpha2_index =[]; 6alpha3_index =[]; 7desplazamientos_calculo =[]; 8coordenadas_calculo =[]; 9singular_FE_sidelength_index =[]; 10 new_entrance_inp_index=[]; 11 end Script 9.4 Step 01 script. step 02 1function [nofelements, nofcrowns, tipo1 , tipo2 , separadores ] = step_02_obtain_parameters (id) 2%% This code finds the crown type data by the naming code 3separadores= strfind (id , _); 4nofelements=str2num(id(1: separadores (1)−2)); 5if str2num(id( separadores (end)+1:end))~=0 6nofcrowns=2; 7else 8nofcrowns=1; 9end 10 tipo1=str2num(id( separadores (2)+1:separadores (3)−1)); 11 tipo2=str2num(id( separadores (4)+1:end)); Script 9.5 Step 02 script. step 03 1function [alpha1_index,alpha2_index,alpha3_index, singular_FE_sidelength_index , file_index ,name_index,ncasos] = step_03_name_archives(nofelements ,nofcrowns,tipo1 , tipo2 , hex_sidelength , id ,name_index, file_index , singular_FE_sidelength_index ,alpha1_index,alpha2_index,alpha3_index, separadores ) 2cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 3 4new_entrance=[]; 5 6%% 1 Corona: 7 8%% tipo_1_tipo_0, tipo_2_tipo_0 , tipo_3_tipo_0 , 9if nofcrowns==1 && (tipo1==1 || tipo1==2 || tipo1==3) 10 alpha1_in=0.05; 11 alpha1_in_plus=0; 9.1 Codes for study 1 129 12 alpha2_in_n=1/cos(pi/nofelements); 13 alpha2_in_plus=0.05; 14 end 15 if nofcrowns==1 && tipo1==4 16 alpha1_in=0.05; 17 alpha1_in_plus=0; 18 alpha2_in_n=1/cos(pi /2/ nofelements); 19 alpha2_in_plus=0; 20 end 21 if nofcrowns==1 && tipo1==5 && (nofelements==4 || nofelements==6 || nofelements==8) 22 alpha1_in=0.05; 23 alpha1_in_plus=0; 24 alpha2_in_n=1/cos(pi /2/ nofelements); 25 alpha2_in_plus=0.0001; 26 end 27 if nofcrowns==1 && tipo1==5 && (nofelements==12 || nofelements==24) 28 alpha1_in=0.05; 29 alpha1_in_plus=0; 30 alpha2_in_n=1/cos(pi /2/ nofelements); 31 alpha2_in_plus=0.0001; 32 end 33 34 %% 2 Coronas: 35 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && (tipo2==1 || tipo2==2 || tipo2==3) 36 alpha1_in=0.05; 37 alpha1_in_plus=0; 38 alpha2_in_n=1/(cos(pi/nofelements)); 39 alpha2_in_plus=0.05; 40 alpha3_in_n=1/(cos(pi/nofelements)); 41 alpha3_in_plus=0.05; 42 end 43 44 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && tipo2==4 45 alpha1_in=0.05; 46 alpha1_in_plus=0; 47 alpha2_in_n=1/(cos(pi/nofelements)); 48 alpha2_in_plus=0; 49 alpha3_in_n=1/(cos(pi /2/ nofelements)); 50 alpha3_in_plus=0; 51 end 52 53 if nofcrowns==2 && (tipo1==1 || tipo1==2 || tipo1==3) && tipo2==5 54 alpha1_in=0.05; 55 alpha1_in_plus=0; 56 alpha2_in_n=1/(cos(pi/nofelements)); 57 alpha2_in_plus=0.05; 58 alpha3_in_n=1/(cos(pi /2/ nofelements)); 59 alpha3_in_plus=0.05; 60 end 61 62 if nofelements==4 && nofcrowns==2 && (tipo1==3) && tipo2==5 63 alpha1_in=0.05; 64 alpha1_in_plus=0; 65 alpha2_in_n=1/(cos(pi/nofelements)); 66 alpha2_in_plus=0; 67 alpha3_in_n=1/(cos(pi /2/ nofelements)); 68 alpha3_in_plus=0.0001; 69 end 70 71 if nofcrowns==2 && (tipo1==4) && (tipo2==1 || tipo2==2 || tipo2==3) 72 alpha1_in=0.05; 73 alpha1_in_plus=0; 74 alpha2_in_n=1/(cos(pi /2/ nofelements)); 75 alpha2_in_plus=0; 76 alpha3_in_n=1/(cos(pi/nofelements)); 77 alpha3_in_plus=0.05; 78 end 79 80 if nofcrowns==2 && (tipo1==4) && (tipo2==4) 81 alpha1_in=0.05; 82 alpha1_in_plus=0; 130 Chapter 9. Appendix: Codes 83 alpha2_in_n=1/(cos(pi /2/ nofelements)); 84 alpha2_in_plus=0; 85 alpha3_in_n=1/(cos(pi /2/ nofelements)); 86 alpha3_in_plus=0; 87 end 88 89 if nofcrowns==2 && (tipo1==4) && (tipo2==5) 90 alpha1_in=0.05; 91 alpha1_in_plus=0; 92 alpha2_in_n=1/(cos(pi /2/ nofelements)); 93 alpha2_in_plus=0; 94 alpha3_in_n=1/(cos(pi /2/ nofelements)); 95 alpha3_in_plus=0; 96 end 97 98 if nofcrowns==2 && (tipo1==5) && (tipo2==1 || tipo2==2 || tipo2==3) && (nofelements==4 || nofelements==6) 99 alpha1_in=0.05; 100 alpha1_in_plus=0; 101 alpha2_in_n=1/(cos(pi /2/ nofelements)); 102 alpha2_in_plus=0.0001; 103 alpha3_in_n=1/(cos(pi/nofelements)); 104 alpha3_in_plus=0; 105 end 106 107 if nofcrowns==2 && (tipo1==5) && (tipo2==1 || tipo2==2 || tipo2==3) && (nofelements==8) 108 alpha1_in=0.05; 109 alpha1_in_plus=0; 110 alpha2_in_n=1/(cos(pi /2/ nofelements)); 111 alpha2_in_plus=0.0001+0.05; 112 alpha3_in_n=1/(cos(pi/nofelements)); 113 alpha3_in_plus=0.05; 114 end 115 116 if nofcrowns==2 && (tipo1==5) && (tipo2==1 || tipo2==2 || tipo2==3) && (nofelements==12) 117 alpha1_in=0.2; 118 alpha1_in_plus=0; 119 alpha2_in_n=1/(cos(pi /2/ nofelements)); 120 alpha2_in_plus=0.0001+0.05; 121 alpha3_in_n=1/(cos(pi/nofelements)); 122 alpha3_in_plus=0.05; 123 end 124 125 if nofcrowns==2 && (tipo1==5) && (tipo2==1 || tipo2==2 || tipo2==3) && (nofelements==24) 126 alpha1_in=0.2; 127 alpha1_in_plus=0; 128 alpha2_in_n=1/(cos(pi /2/ nofelements)); 129 alpha2_in_plus=0.0001+0.05; 130 alpha3_in_n=1/(cos(pi/nofelements)); 131 alpha3_in_plus=0; 132 end 133 134 if nofcrowns==2 && (tipo1==5) && (tipo2==4) && nofelements~=12 && nofelements~=24, 135 alpha1_in=0.05; 136 alpha1_in_plus=0; 137 alpha2_in_n=1/(cos(pi /2/ nofelements)); 138 alpha2_in_plus=0.0001+0.05; 139 alpha3_in_n=1/(cos(pi /2/ nofelements)); 140 alpha3_in_plus=0.05; 141 end 142 143 if nofcrowns==2 && (tipo1==5) && (tipo2==4) && (nofelements==12) 144 alpha1_in=0.2; 145 alpha1_in_plus=0; 146 alpha2_in_n=1/(cos(pi /2/ nofelements)); 147 alpha2_in_plus=0.0001; 148 alpha3_in_n=1/(cos(pi /2/ nofelements)); 149 alpha3_in_plus=0.05; 150 end 151 152 if nofcrowns==2 && (tipo1==5) && (tipo2==4) && (nofelements==24) 153 alpha1_in=0.2; 9.1 Codes for study 1 131 154 alpha1_in_plus=0; 155 alpha2_in_n=1/(cos(pi /2/ nofelements)); 156 alpha2_in_plus=0.0001+0.05; 157 alpha3_in_n=1/(cos(pi /2/ nofelements)); 158 alpha3_in_plus=0.05; 159 end 160 161 if nofcrowns==2 && (tipo1==5) && (tipo2==5) && nofelements~=12 && nofelements~=24, 162 alpha1_in=0.05; 163 alpha1_in_plus=0; 164 alpha2_in_n=1/(cos(pi /2/ nofelements)); 165 alpha2_in_plus=0.0001+0.05; 166 alpha3_in_n=1/(cos(pi /2/ nofelements)); 167 alpha3_in_plus=0.05; 168 end 169 170 if nofcrowns==2 && (tipo1==5) && (tipo2==5) && (nofelements==12) 171 alpha1_in=0.2; 172 alpha1_in_plus=0; 173 alpha2_in_n=1/(cos(pi /2/ nofelements)); 174 alpha2_in_plus=0.0001+0.05; 175 alpha3_in_n=1/(cos(pi /2/ nofelements)); 176 alpha3_in_plus=0.05; 177 end 178 179 if nofcrowns==2 && (tipo1==5) && (tipo2==5) && (nofelements==24) 180 alpha1_in=0.2; 181 alpha1_in_plus=0; 182 alpha2_in_n=1/(cos(pi /2/ nofelements)); 183 alpha2_in_plus=0.0001+0.05; 184 alpha3_in_n=1/(cos(pi /2/ nofelements)); 185 alpha3_in_plus=0.05; 186 end 187 188 %% One−crown 189 if nofcrowns==1, 190 for alpha1=alpha1_in+alpha1_in_plus :0.05:0.85 191 for alpha2=alpha1∗alpha2_in_n+alpha2_in_plus :0.05:0.85 192 singular_FE_sidelength =alpha1∗hex_sidelength; 193 alpha1_s=num2str((alpha1) ,4) ; 194 alpha1_s=strrep (alpha1_s, '.',' ' ) ; 195 alpha2_s=num2str((alpha2), '%2.15f'); 196 alpha2_s=strrep (alpha2_s, '.',' ' ) ; 197 alpha2_s=alpha2_s(1:4) ; 198 Files=dir ([ '∗' ,'_',alpha1_s, '_',alpha2_s(1:end), ' ∗.inp ']) ; 199 if isempty( Files ) 200 alpha2_s=strrep (num2str(alpha2+0.0001,'%2.15f') , '.',' ' ); 201 alpha2_s=alpha2_s(1:4) ; 202 Files =[]; 203 Files=dir ([ '∗' ,'_',alpha1_s, '_',alpha2_s(1:end), ' ∗.inp ']) ; 204 if isempty( Files ) 205 alpha2_s=strrep (num2str(alpha2−0.0001,'%2.15f') , '.',' ' ); 206 alpha2_s=alpha2_s(1:4) ; 207 end 208 end 209 new_entrance=[id alpha1_s '_'alpha2_s ]; 210 new_entrance_inp=[id '_'alpha1_s '_'alpha2_s '.inp ']; 211 name_index=strvcat(name_index, new_entrance); 212 file_index = strvcat ( file_index ,new_entrance_inp); 213 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; singular_FE_sidelength ]; 214 alpha1_index=[alpha1_index,alpha1 ]; 215 alpha2_index=[alpha2_index,alpha2 ]; 216 end 217 end 218 end 219 220 %% Two−crowns 221 if nofcrowns==2 222 for alpha1=alpha1_in+alpha1_in_plus :0.05:0.85 223 for alpha2=alpha1∗alpha2_in_n+alpha2_in_plus :0.05:0.85 224 for alpha3=alpha2∗alpha3_in_n+alpha3_in_plus :0.05:0.85 132 Chapter 9. Appendix: Codes 225 % Singular element side−length 226 singular_FE_sidelength =alpha1*hex_sidelength; 227 228 % alpha1 229 alpha1_s=num2str((alpha1) ,4) ; 230 alpha1_s=strrep (alpha1_s, ., ); 231 232 % alpha2 233 alpha2_s=num2str((alpha2), %2.15f ); 234 alpha2_s=strrep (alpha2_s, ., ); 235 alpha2_s=alpha2_s(1:4) ; 236 237 % alpha3 238 alpha3_s=num2str((alpha3), %2.15f ); 239 alpha3_s=strrep (alpha3_s, ., ); 240 if length(alpha3_s)>3, alpha3_s=alpha3_s(1:4) ; else , alpha3_s=alpha3_s(1:end); end 241 242 % = alpha2, = alpha3 243 Files=dir ([ *,_,alpha1_s, _,alpha2_s(1:end), _,alpha3_s(1:end), *.inp ]) ; 244 if isempty( Files ), 245 % = alpha2, + alpha3. 246 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 247 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 248 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 249 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 250 alpha3_saux=alpha3_saux(1:4); 251 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 252 if ~isempty( Files ) , 253 alpha2_s=alpha2_saux; 254 alpha3_s=alpha3_saux; 255 else 256 % = alpha2, −alpha3. 257 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 258 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 259 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 260 alpha3_saux=alpha3_saux(1:4); 261 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 262 if ~isempty( Files ) , 263 alpha2_s=alpha2_saux; 264 alpha3_s=alpha3_saux; 265 else 266 % + alpha2, = alpha3 267 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 268 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 269 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 270 alpha2_saux=alpha2_saux(1:4); 271 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 272 if ~isempty( Files ) , 273 alpha2_s=alpha2_saux; 274 alpha3_s=alpha3_saux; 275 else 276 % + alpha2, + alpha3 277 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 278 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 279 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 280 alpha2_saux=alpha2_saux(1:4); 281 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 282 alpha3_saux=alpha3_saux(1:4); 283 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 284 if ~isempty( Files ) , 285 alpha2_s=alpha2_saux; 286 alpha3_s=alpha3_saux; 287 else 288 % + alpha2, − alpha3 289 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 290 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 291 alpha2_saux=strrep (num2str(alpha2_aux+0.001, %2.15f ) , ., ); 292 alpha2_saux=alpha2_saux(1:4); 293 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 294 alpha3_saux=alpha3_saux(1:4); 295 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 9.1 Codes for study 1 133 296 if ~isempty( Files ) , 297 alpha2_s=alpha2_saux; 298 alpha3_s=alpha3_saux; 299 else 300 % − alpha2, = alpha3 301 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 302 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 303 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 304 alpha2_saux=alpha2_saux(1:4); 305 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 306 if ~isempty( Files ) , 307 alpha2_s=alpha2_saux; 308 alpha3_s=alpha3_saux; 309 else 310 % − alpha2, + alpha3 311 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 312 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 313 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 314 alpha2_saux=alpha2_saux(1:4); 315 alpha3_saux=strrep (num2str(alpha3_aux+0.001, %2.15f ) , ., ); 316 alpha3_saux=alpha3_saux(1:4); 317 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 318 if ~isempty( Files ) , 319 alpha2_s=alpha2_saux; 320 alpha3_s=alpha3_saux; 321 else 322 % − alpha2, − alpha3 323 alpha2_aux=str2num([alpha2_s(1) , .,alpha2_s(2:end)]) ; 324 alpha3_aux=str2num([alpha3_s(1) , .,alpha3_s(2:end)]) ; 325 alpha2_saux=strrep (num2str(alpha2_aux−0.001, %2.15f ) , ., ); 326 alpha2_saux=alpha2_saux(1:4); 327 alpha3_saux=strrep (num2str(alpha3_aux−0.001, %2.15f ) , ., ); 328 alpha3_saux=alpha3_saux(1:4); 329 Files=dir ([ *,_,alpha1_s, _,alpha2_saux(1:end), _,alpha3_saux(1:end), *.inp ]) ; 330 if ~isempty( Files ) , 331 alpha2_s=alpha2_saux; 332 alpha3_s=alpha3_saux; 333 end 334 end 335 end 336 end 337 end 338 end 339 end 340 end 341 end 342 if isempty( Files ) 343 new_entrance_inp=[]; 344 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; []]; 345 name_index=strvcat(name_index, new_entrance); 346 alpha1_index=[alpha1_index ,[]]; 347 alpha2_index=[alpha2_index ,[]]; 348 alpha3_index=[alpha3_index ,[]]; 349 else 350 new_entrance_inp=[id _alpha1_s _alpha2_s _alpha3_s .inp ]; 351 name_index=strvcat(name_index, new_entrance); 352 file_index = strvcat ( file_index ,new_entrance_inp); 353 singular_FE_sidelength_index =[ singular_FE_sidelength_index ; singular_FE_sidelength ]; 354 alpha1_index=[alpha1_index,alpha1 ]; 355 alpha2_index=[alpha2_index,alpha2 ]; 356 alpha3_index=[alpha3_index,alpha3 ]; 357 end 358 end 359 end 360 end 361 end 362 363 ncasos=length( file_index ); Script 9.6 Step 03 Script. 134 Chapter 9. Appendix: Codes 9.1.3 step 04 1function [coordenadas,coordenadasN5,coordenadas5, triangulos , triangulos5 ] = step_04_mesher(inp_malla,nofelements, sidelength ,alpha,typeelement, vct); 2xcentro=vct(1) ; ycentro=vct(2) ; zcentro=vct(3) ; 3 4%% Matriz de coordenadas 5begin_reading_coordenadas= strfind (inp_malla, *Node ); % Leo desde "*Node". 6end_reading_coordenadas= strfind (inp_malla, *Element ) ; % Hasta "*Element". 7coordenadas_aux=[inp_malla(begin_reading_coordenadas+8:end_reading_coordenadas−2)]; % Saco las coordenadas del documento. El +5 es para leer el inicio de los números, y no el encabezado. El −2 es para no coger los caracteres del encabezado del siguiente punto. 8coordenadas_aux= strsplit (coordenadas_aux, \n ); % Corto el texto por l íneas según el caracter " salto de l ínea". 9coordenadas_aux=coordenadas_aux ; % Ordeno el texto en vertical . 10 for i=1:length(coordenadas_aux) 11 coordenadas(i ,:) =str2num(cell2mat(coordenadas_aux(i))); % Convierto cada fila en matriz de texto , y cada matriz de texto ( cell ) en matriz numérica. 12 end 13 coordenadas (:,1) =[]; % Elimino la columna del número de identificaci ón del nodo. 14 15 %% Matriz de conectividad de elementos 16 17 begin_reading_triangulos = strfind (inp_malla, *Element, ) ; % Se encuentra el comienzo de la matriz de tri ángulos 18 end_reading_triangulos = strfind (inp_malla, *End Part ) ; % Se encuentra el final de la matriz de tri ángulos 19 triangulos_aux =[inp_malla( begin_reading_triangulos +21: end_reading_triangulos −2)]; % Se introducen los valores en una variable 20 triangulos_aux = strsplit ( triangulos_aux , \n ); % Se divide las l íneas de la matriz por el caracter de salto de l ínea 21 triangulos_aux =triangulos_aux ; % Se pone en el formato correcto 22 for i=1:length( triangulos_aux ) 23 triangulos (i ,:) =str2num(cell2mat( triangulos_aux ( i))); % Se convierte la l ínea de la variable tri ángulos en nú meros 24 end 25 triangulos (:,1) =[]; % Se elimina la primera columna, que es innecesaria 26 27 %% Matriz de coordenadas de los elementos singulares 28 29 N=2*nofelements; % El número de lados del polígono regular cuya mitad superior ser ían los elementos singulares . 30 n=(0:nofelements); 31 r=sidelength *1/2; 32 if typeelement==1 % Si existen elementos singulares 33 x(n+1)=r*cos(2*pi*n/N); 34 y(n+1)=r*sin(2*pi*n/N); % Calculo las coordenadas x, y del polígono regular . 35 z(n+1)=zcentro; % Introducir la z pertinente . 36 coordenadas_singulares =[xcentro , ycentro , zcentro ; x , y , z ]; 37 elseif typeelement==0 % Si no existen los elementos singulares 38 coordenadas_singulares =[xcentro , ycentro , zcentro ]; 39 end 40 coordenadas_singulares (:,3) =[]; % Se elimina la tercera columna, innecesaria 41 42 %% Matriz de conectividad de los elementos singulares 43 44 r=sidelength ; 45 x(n+1)=r*cos(2*pi*n/N); 46 y(n+1)=r*sin(2*pi*n/N); 47 z(n+1)=0; % Introducir la z pertinente . 48 for n=1:nofelements+1 49 closestpoint (n)=dsearchn(coordenadas,[x(n) y(n)]) ; % Encuentro el punto del listado más próximo al obtenido de nuestra fórmula 50 end 51 if typeelement==1 % Se crea la matriz de conectividad de los elementos singulares . Los puntos están en el sentido de las agujas del reloj : el primer nodo es el central , el segundo nodo es el siguiente al central , el tercero es el punto que ya exist ía previamente en la malla (el final del elemento singular ) , el cuarto el punto siguiente al anterior (que también preexiste en la malla), y el último el siguiente punto intermedio . 52 for i=1:nofelements 53 triangulos5 (i ,:) =[size (coordenadas,1)+1 size (coordenadas,1)+1+i closestpoint (i) closestpoint (i+1) size ( coordenadas,1)+2+i]; 54 end 55 elseif typeelement==0 56 for i=1:nofelements % Aquí no hay puntos intermedios , solo tres 57 triangulos5 (i ,:) =[size (coordenadas,1)+1 closestpoint (i) closestpoint (i+1)]; 9.1 Codes for study 1 135 58 end 59 triangulos =[ triangulos ; triangulos5 ]; 60 triangulos5 =[]; 61 end 62 %% Terminamos de ensamblar la matriz de coordenadas 63 coordenadas5=coordenadas_singulares ; 64 coordenadasN5=coordenadas; 65 coordenadas=[coordenadas; coordenadas_singulares ]; Script 9.7 Step 03 Script. step 13 1function [ ajuste , ajuste_2 , slope_index ,max_2nd_crown_change,casos_buenos,nota_a,nota_b, criterio_1 , criterio_2 , criterio_3a , criterio_3b , criterio_4 ,criterio_1_num, criterio_1_den , criterio_5 , Coeff_vectora , Coeff_vectorb , mean_slopes_index, criterio_2_mejorado ] = step_13_numerical_regressions (alpha1_index,alpha2_index,alpha3_index , desplazamientos_calculo , coordenadas_calculo , tipo1,tipo2 , id, separadores ,coordenadas, rcts , t ) 2cd([ C:\Users\Usuario\Desktop\TFG (Avances y Bibliografia )\Macros\ , id (1: separadores (1)−1), \, id ]) ; 3% clc 4 5%% X: alpha_2, Y: alpha_3, Z: u_ct; X coordinate : alpha_2; Y coordinate : alpha_3; Z coordinate : u_ct 6 7%% Variables initialization 8leyenda =[]; ajuste =[]; 9x_pr=[]; y_pr=[]; z_pr =[]; 10 casos_buenos=[]; slope_index =[]; level_curve =[]; point_vector =[]; max_2nd_crown_change=[]; mean_slopes_index=[]; slope_index_aux=[]; criterio_2_mejorado =[]; 11 change=[1]; change2=[1]; 12 13 %% This section fills sparce matrices (FEM results). 14 desplazamientos_calculo =full ( desplazamientos_calculo ); 15 coordenadas_calculo=full ( coordenadas_calculo); 16 17 %% This section finds where the indices change, that is , all the cases with the same alpha_1 value. 18 19 20 for s=1:length (alpha1_index)−1 21 if alpha1_index(s)−alpha1_index(s+1)~=0 % Each value that does not coincide with the following 22 change=[change, s+1]; % its position is added as the first unequal values . 23 end 24 end 25 26 %% For the one−crown cases: 27 if tipo2==0 28 29 %% Finding an optimal point . We must find the right proportion between alpha_1 and alpha_2 that results in an u_ct value of 2.25596. To do that , we plot u_ct vs alpha_2 for each alpha_1 ( there will be several curves) . If we also plot u_ct=2.25596 we can find a cut−off point between curves (since they are defined by points , we will be using a spline to find the cut−off point) . For each alpha_1 value we find the alpha_2 and alpha_3 values that gives the right u_ct. A simple numerical correlation is presupposed. 30 % For one−crown cases, we can scatter alpha_1 vs optimal alpha_2 values . Then, a ployfit regression is made. 31 32 for s2=1:length (change)−1 % For each curve ( there are as many curves as changes in vector alpha_1) taken as an entry of the vector change 33 34 %% The crack tip displacements are plotted 35 alpha1_legend=num2str(alpha1_index(change(s2))); 36 % First figure : u_ct vs alpha_2 for each alpha_1 37 figure (1) 38 xlabel (\alpha_2 ); ylabel (u_1 ); 39 plot (alpha2_index(change(s2):change(s2+1)−1), desplazamientos_calculo (change(s2):change(s2+1)−1,1), linewidth ,2) 40 if s2~=length(change)−1 41 leyenda= strvcat (leyenda ,[ \alpha_1= ,alpha1_legend]) ; 42 end 43 hold on 44 45 %% This section calculates the intersection between u_ct and alpha_1. A spline is generated with the previously calculated points (the FEM results).