In a linear elastic analysis, it is well known that, at the tip of a crack the stress field becomes infinite and thus, is singular. The strength of this singularity is measured by the SIF that is thus defined at the crack tip. The presence of the stress singularity in the numerical model raises considerable numerical difficulties, by virtue of the need of simultaneously representing the singular and the finite stresses in the numerical model. Instead of representing the stress singularity in the numerical model, Oliveira and Portela [26] used an elegant strategy that subtracts the singularity from the elastic field which is known as the singularity subtraction technique (SST). Hence, the SST performs a regularization of the stress field, which introduces the SIF as primary unknowns of the numerical method used in the analysis. These two features, which are the analysis of the regularized stress field and the direct computation of SIF, make very efficient the SST solution strategy. The paper considers the SST, a very efficient and accurate technique for solving twodimensional problems of linear elastic fracture mechanics, as reported by Oliveira et al. [27], implemented in the ILMF mesh-free model of numerical analysis.

Mesh-free numerical methods eliminate the mesh of the discretization, an intrinsic feature of the finite element and finite difference numerical methods of the first-generation in computational mechanics. On the other hand, the development of the boundary element method, as a second-generation numerical method, was motivated by the reduction of the analysis dependency on the mesh discretization, I. Introduction used only on the boundary of the domain. Mesh-free methods are third-generation numerical methods which consider only a nodal discretization and completely overcome the difficulties posed by the mesh of the first and second-generation numerical methods in computational mechanics. This paper considers a domain mesh-free method of analysis, with the MLS approximation of the elastic field, coupled with a multi-objective optimization process that automatically generates optimal nodal arrangements of the mesh-free discretization, to compute the SIF of two-dimensional linear elastic fracture mechanics problems.

Thorough reviews of mesh-free methods and their applications in science and engineering were recently presented by Chen et al. [6] and Huerta et al. [17]. The most popular of local mesh-free methods is the MLPG method, presented by Atluri and Zhu [2] to Atluri and Shen [1], which implements the MLS approximation. Other local mesh-free methods of reference are the LPIM method, see Liu and Gu [21] and the LRPIM method, see Liu et al. [22]. The ILMF linearly integrated local mesh-free method, presented by Oliveira et al. [27], performs a linear reduced integration, which leads to an increase of the solution accuracy with high efficiency.

Until now, the discretization process of local mesh-free methods has been heuristically implemented, which requires an expensive and time consuming calibration of the nodal arrangements or parameters of the discretization that refer to the size of the compact supports and the size of the integration domain of each node. This is a huge drawback since the definition of these discretization parameters is not unique and therefore cannot be easily implemented into an automatic process.

Some researchers tried to overcome the drawback of heuristically defined meshfree discretization parameters, as is the case of Baradaran and Mahmoodabadi [4], Bagheri et al. [3] and Ebrahimnejad et al. [12]. The successful attempts of these authors required an analytical solution to be performed and therefore their modeling process is not efficient. Recently, Santana et al. [31] and Oliveira and Portela [26], presented a strategy that performs the optimization of the size of compact supports and the size of the local integration domain of given mesh-free nodal distributions.

Thus, there is room for the alternative modeling strategy of this paper, that is the automatic generation of optimal mesh-free parameters and nodal discretization, through an optimization process that completely overcomes the issues of the heuristic process of discretization. As a consequence, the modeling strategy of this paper ensures robustness, accuracy and efficiency of the analysis, features required to be able to make reliable statements in the high fidelity modeling of engineering applications.

The use of optimization has been applied in many different areas, such as elastostatics, see Denk et al. [10], Proos et al. [29] and Zolfagharian et al. [37], heat conduction, see Dede [8], Denk et al. [11], Gersborg-Hansen et al. [14] and Kim et al. [20], fluid mechanics, see Dede et al. [9], electrostatics, see Gupta et al. [15], or structural dynamics, see Kim et al. [20] and Proos et al. [29].

The field of optimization is expansive, and the choice of a suitable algorithm is highly problem dependent, as reported by Zingg et al. [36]. The No free lunch theorems for optimization, presented by Wolpert and Macready [35], suggests that different algorithms are better than others for particular classes of problems. The multi-objective optimization of mesh-free numerical models deals with two main difficulties. The first one concerns the number of optimal solutions, generated by competing goals, instead of a single optimal solution. The second difficulty regards the large and complex search space that cannot be dealt with classic optimization methods. Consequently, to overcome these difficulties, non-gradient methods of optimization must be used, instead of classic methods. Evolutionary algorithms are non-gradient methods, quite robust in locating the global optimum, that do not require continuity or predictability over the design space. This paper considers the use of evolutionary genetic algorithms (GA), for the multi-objective optimization of nodal arrangements of the mesh-free discretization. GA perform a search and optimization procedure motivated by the principles of natural genetics and natural selection, originally proposed by Holland [16]. They are a robust and flexible approach that can be applied to a wide range of optimization problems, as as reported for instance by Kelner and Leonard [19], McCall [23] and Ebrahimnejad et al. [12].

The paper organization is as follows. The modeling of the structural body and the local mesh-free method is presented in Section 2 that is followed by the implementation of the SST in the mesh-free formulation, presented in Section 3. Section 4 presents the multi-objective optimization implementation and algorithm formulation. Numerical results, obtained for benchmark problems, in order to illustrate the accuracy, efficiency and robustness of the strategies adopted in this work, are presented in Section 5. Finally, the concluding remarks are presented in Section 6.

The local mesh-free numerical analysis of the structural body is carried out by the ILMF model presented by Oliveira et al. [27]. It is defined in a body with domain ? and boundary Î?" = Î?" u ? Î?" t , with constrained displacements u prescribed on the kinematic boundary Î?" u and loaded by an external system of distributed surface and body forces, with densities represented respectively by t, applied on the static boundary Î?" t , and b, applied in ?, as Figure 1 schematically represents. Assign to Mesh-free discretization of a body with domain ? and boundary Î?" = Î?" u ? Î?" t ; ? P , ? Q and ? R are local domains assigned to reference nodes P , Q and R; The work theorem is used to formulate the ILMF model. The mechanical equilibrium of the local domain ? Q can be defined through the rigid-body kinematic formulation of the work theorem, as presented by Oliveira et al. [27], which is II. Mesh-free Modeling of the Structural Body written, in the case of no body forces, as
? Q has boundary Î?" Q = Î?" Qi ? Î?" Qt ? Î?" Qu , in which Î?" Qi is the interior local boundary and Î?" Qt ? Î?" t and Î?" Qu ? Î?" u . point Q an arbitrary local domain ? Q , such that Q ? ? Q ? ? ? Î?", with boundary Î?" Q = Î?" Qi ? Î?" Qt ? Î?" Qu ,Î?" Q ?Î?" Qt t dÎ?" = ? Î?" Qt t dÎ?" (1)
and describes the equilibrium of boundary tractions in ? Q . This equation, used to generate the stiffness matrix of each node of a mesh-free discretization, is integrated by Gauss quadrature. Finally, in order to allow for a unique solution of the elastic field, displacement boundary conditions must be enforced, on the kinematic boundary Î?" u , as
u = u. (2)
Since ILMF is a local model, each node of the discretization has assigned its local integration domain, as schematically represented in Figure 1, which has rectangular or circular shape, as Figure 2 schematically represents. Whenever a linear q (a) Rectangular q (b) Circular Schematic representation of local integration domains, with 1 integration point per boundary, or quadrant, of the local domain, for the computation of local equilibrium equations. variation of tractions is defined, along each segment of the boundary of the local domain, equilibrium equations (1) can be exactly evaluated with 1 integration point, centered on each segment of the local boundary, as represented in Figure 2, which leads to a point-wise discrete form of mechanical equilibrium, represented by
L i n i n i j=1 t x j = ? L t n t nt k=1 t x k ,(3)
where the number of integration points, or segments, defined on, respectively the interior boundary Î?" Qi = Î?" Q ? Î?" Qt ? Î?" Qu , of length L i , and the static boundary Î?" Qt , of length L t , are denoted, respectively by n i and n t . The MLS approximation of variables is used with the ILMF model. Therefore, traction components are evaluated in terms of the unknown nodal parameters û and thus, equations (1) lead to the system of algebraic equations of the order 2 × 2n (n is the number of nodes of the influence domain of the reference node Q), given by
Î?" Q ?Î?" Qt nDB dÎ?" û = ? Î?" Qt t dÎ?"(4)
represented by in which K Q is the stiffness matrix
K Q û = F Q ,(5)K Q = Î?" Q ?Î?" Qt nDB dÎ?" (6)
and F Q is the vector of forces
F Q = ? Î?" Qt t dÎ?".(7)
For a nodal arrangement with N nodes, in which M are interior and static-boundary nodes, the assembly of equations ( 5) leads to the system of equations of the order
2M × 2N K û = F.(8)
The remaining algebraic equations are generated from the N ? M kinematicboundary nodes, through the direct interpolation of the boundary condition (2) prescribed as
u k = ? k û = u k ,(9)
with k = 1, 2, where u k is the constrained displacement component. Equations ( 9) are directly assembled into the global system (8).

For each node of a local mesh-free discretization there are two key parameters, respectively the size r ?s of the compact support ? s and the size r ?q of the local integration domain ? q that strongly affect the performance of the solution. For a generic node i, these parameters are defined through arbitrary constants, ? s and ? q , respectively as r ?s = ? s c i (10) and
r ?q = ? q c i ,(11)
in which c i is the distance of the node i to the nearest neighboring node. Equations (10) and (11) show that the accuracy of a mesh free numerical application can be controlled through a proper specification of the discretization parameters ? s and ? q .

It is important to enhance the different roles that these parameters play, in any numerical application. The size of the influence domain of a point, determined by the size of the compact support of each node, completely defines the number of nodes used to build MLS shape functions of that point. Consequently, the parameter ? s , sometimes referred to as MLS discretization parameter, is primarily linked to the accuracy of the numerical application. On the other hand, since the integration domain of each node is used to compute the respective nodal stiffness matrix, it must be entirely defined within the domain of the body, without intersecting the respective boundary. Consequently, the parameter ? q , sometimes referred to as the local domain parameter, is linked primarily to the efficiency of the application.

Until now, these parameters have been heuristically defined with values that depend on the pattern of the nodal distribution of each mesh-free application. When changes are made to the nodal distribution of a problem, these parameters also change, becoming an interactive process. In general, they have been considered in the range, respectively of ? s > 1.0 and ? q < 1.0, as reported by Oliveira et al. [27]. In this paper, the appropriate values of ? s , ? q and the nodal distribution of nodes are obtained automatically, through a multi-objective optimization process.

To overcome difficulties raised by the presence of unbounded stress in the numerical method, an alternative strategy which considers the subtraction of singularities from the original elastic field is used. The leading feature of this formulation is the regularization the elastic field, through the subtraction of the crack tip singularity from the original elastic field, which introduces the SIF as additional primary unknowns of the regularized numerical model.

In the linear elastic fracture mechanics, the stress field is singular at a crack tip and therefore, it is convenient to modify the original problem before its solution by the ILMF numerical model. As the linear behavior allows the principle of superposition, the elastic field can be decomposed into a regular (R) and a singular (S) component as
? ij = (? ij ? ? S ij ) + ? S ij = ? R ij + ? S ij(12)
and
u i = (u i ? u S i ) + u S i = u R i + u S i ,(13)
where ? R ij = ? ij ?? S ij and u R i = u i ?u S i represent the regular parts, respectively of the stress and displacement of the initial problem; ? S ij and u S i denote, respectively the stress and displacement of a particular solution, of the initial problem, which represent the singular field. When suitable functions are used for the particular singular field, equations ( 12) and ( 13) regularize the initial problem, since the stress ? R ij become non-singular. With this regularization, the analysis of the initial problem can be performed with the regular elastic field only, since the components ? S ij and u S i automatically satisfy the field equations identically, because they represent a particular solution of the initial problem. Therefore, the elasticity equations are written as
L T ? R = 0 (14) ? R = L u R (15) ? R = D ? R(16)
in domain ?, with boundary conditions
u R = u ? u S on Î?" u (17) and t R = t ? t S on Î?" t .(18)
Note that the boundary conditions ( 17) and ( 18) include additional terms, respectively u S and t S , components of a singular particular solution of the initial problem.

Components ? S ij and u S i of the particular solution, used in equations ( 12) and ( 13), represent the singular field around the crack tip, which can be defined through the first term of the William's [34] eigen-expansion, derived for a semi-infinite edge crack. The stress components are
III. The Singularity Subtraction Technique -SST a) Regularized Elastic Field b) William's Singular Solution ? S 11 = K I ? 2?r cos ? 2 1 ? sin ? 2 sin 3? 2 + K II ? 2?r sin ? 2 2 + cos ? 2 cos 3? 2 ,(19)? S 22 = K I ? 2?r cos ? 2 1 + sin ? 2 sin 3? 2 ? K II ? 2?r sin ? 2 cos ? 2 cos 3? 2(20)
and
? S 12 = K I ? 2?r cos ? 2 sin ? 2 cos 3? 2 + K II ? 2?r cos ? 2 1 ? sin ? 2 sin 3? 2(21)
and the displacements are
u S 1 = K I 4µ r 2? (2? ? 1) cos ? 2 ? cos 3? 2 + + K II 4µ r 2? (2? + 3) sin ? 2 + sin 3? 2(22)
and
u S 2 = K I 4µ r 2? (2? + 1) sin ? 2 ? sin 3? 2 + + K II 4µ r 2? (2? ? 3) cos ? 2 + cos 3? 2 ,(23)
where K I and K II represent the SIF, respectively of the opening and sliding modes; the constant ? = 3 ? 4? is defined for plain strain and ? = (3 ? ?)/(1 + ?) for plain stress, in which ? is Poisson's ratio; the constant µ is the shear modulus. A polar coordinate reference system (r, ?), centered at the crack tip, is defined such that ? = 0 is the crack axis, ahead of the crack tip. Note that the order r ?1/2 of the stress field becomes singular when r tends to zero. Caicedo and Portela [5] demonstrated that the first term of the William's eigen-expansion, derived for an edge crack, can also be used to represent the elastic field around the crack-tip, for the case of internal piecewise-flat multi-cracked finite plates, under mixed-mode deformation. At a boundary point, the singular stress components, of equations ( 19) to (21), are used in the definition of traction components as
t S = t S 1 t S 2 = ? S 11 ? S 21 ? S 12 ? S 22 n 1 n 2 = g 11 g 12 g 21 g 22 K I K II = g k,(24)
where n i refers to the i-th component of the unit normal to the boundary, outwardly directed; functions g ij = g ij (r ?1/2 , ?) were introduced for a simple notation of equations ( 19) to (21) and the vector k contains the SIF components.

The displacement field, of equations ( 22) and (23), can be similarly defined in a vector form as
u S = u S 1 u S 2 = f 11 f 12 f 21 f 22 K I K II = f k,(25)
where functions f ij = f ij (r 1/2 , ?) are a simple notation of equations ( 22) and (23). An approximate solution of the regularized problem, equations ( 14) to ( 16) with boundary conditions ( 17) and ( 18), obtained with the ILMF numerical model, is now considered. The equilibrium equations (1), of the domain ? Q associated with the node
Q ? ? Q ? Î?" Q , are now rewritten as Î?" Q ?Î?" Qt t R dÎ?" = ? Î?" Qt t ? t S dÎ?",(26)
in which the static boundary conditions (18), of the regularized problem, are considered. For a linear reduced integration, along each boundary segment of the local domain, equation ( 26) simply leads to
L i n i n i j=1 t R x j = ? L t n t nt k=1 t x k + Î?" Qt t S dÎ?",(27)
in which n i and n t denote the total number of integration points, or boundary segments, defined on, respectively the interior local boundary  27) is done with the MLS approximation, see Oliveira et al. [25], in terms of the unknown nodal parameters ûR , which leads to the system of two linear algebraic equations
Î?" Qi = Î?" Q ? Î?" Qt ? Î?" Qu ,L i n i n i j=1 n x j DB x j ûR = ? L t n t nt k=1 t x k + Î?" Qt g dÎ?" k(28)
that can be written as
K Q ûR + G Q k = F Q ,(29)
in which the stiffness matrix K Q , of the order 2 × 2n (n is the number of nodes included in the influence domain of the node Q) is given by
K Q = L i n i n i j=1 n x j DB x j ,(30)
matrix G Q , of the order 2 × 2, computed from equations (24), is given by
G Q = ? Î?" Qt g dÎ?"(31)
and F Q is the force vector given by
F Q = ? L t n t nt k=1 t x k . (32)
Note that, in the case of an interior node, matrix G Q and vector F Q are null. For a problem with N nodes, the assembly of equations equations ( 29) for all M interior and static-boundary nodes generates the global system of 2M ×(2N +2) equations The N ? M kinematic-boundary nodes, are used to generate the remaining equations of the discretization, implementing the kinematic boundary conditions of the regularized problem, equations (17). Thus, for a kinematic-boundary node, the boundary conditions of the regularized problem are enforced by a direct interpolation method as
K ûR + G k = F.(33u R k = ? k ûR = u k ? u S k = u k ? f k k,(34)
with k = 1, 2, where u k denotes the specified displacement component and u S k = f k k is the displacement component of the singular solution, obtained from equations (25). For the sake of simplicity, equations (34) are written in the same form of equations ( 29), for a point Q, as
K Q k ûR + G Q k k = F Q k ,(35)
in which 35) are assembled into the global system of equations ( 33) which, after this operation, is written as
K Q k = ? k , while G Q k = f k and F Q k = u k . Local equations (K G ûR k = F ,(36)
in which K is a matrix of the order 2N × 2N , G is a matrix of the order 2N × 2 and F is a vector of the order 2N ; the unknowns are the vector ûR , of the order 2N , and the vector k of the order 2. Note that this global system of equations introduce the SIF K I and K II , in the vector k, as additional unknowns of the numerical method. Therefore, to have a well-posed problem, with a unique solution, it is necessary to specify additional constraint equations, one for each mode of deformation considered in the analysis. These additional constraint equations can be specified in two additional bottom rows in the system of equations (36).

The required additional constraints enforce the singularity cancellation in the regularized problem and can be implemented by the cancellation of the regular regular stress components, as
? R ij = 0 ? ? ij = ? S ij(37)
which ensure that, at the crack tip, the initial problem is singular. In order to be effective, the additional constraints must be defined in terms of the unknown regularized nodal parameters of ûR . Conditions (37) can be redefined, in terms of the respective traction components at the crack tip, as
t R j = ? R ij n i = 0 ? t j = t S j ,(38)
where n i denotes the unit normal components of the crack faces. After the MLS approximation, conditions (38), defined at the crack tip x tip , are written as
t R x tip = n x tip D B x tip ûR = 0,(39)
or
C ûR = 0,(40)
in which matrix C = n x tip D B x tip and can now be included in the global system of equations (36), leading to the final system of equations of the order
(2N + 2) × (2N + 2) K G C 0 ûR k = F 0 ,(41)

# Additional Constraints d)

which represents a generalized saddle point problem that can be solved, since the stiffness matrix K, of the ILMF local mesh free model is always non singular, with very low condition numbers, as reported by Oliveira and Portela [26].

The optimization literature contains the basic concepts and terminology required to carry out the optimization process presented in this work, here formally represented by Sawaragi et al. [32], Hwang and Masud [18], Ringuest [30] and Steuer [33].

Multi-objective optimization of the mesh-free model is carried out through an automated procedure that modifies the design or decision variables which are the mesh-free discretization parameters and the nodal distribution. Hence, the optimization process incrementally updates the design variables, carries out a meshfree numerical analysis of the updated model and scans the results of each increment to check if an optimized solution has been reached. In this process, the objective functions define the goal of the optimization, while constraints keep within bounds the value of a design response. The goal of the optimization aims to minimize the objective functions by finding feasible solutions, which are arrangements of mesh-free discretization satisfying the constraints of the problem. It is important to note that the optimizer never deals with solution errors of the generated arrangements of the mesh-free model. Genetic Algorithm (GA) belongs to a class of evolutionary algorithms, defined as a non-derivative global search heuristic, motivated by the principles of natural genetics and natural selection, presented by Holland [16]. GA is an optimization technique that can be applied to a wide range of problems, as seen in Kelner and Leonard [19] and McCall [23], and can also be applied to mesh-free methods, as seen in Bagheri et al. [3] and Ebrahimnejad et al. [12].

The GA keep a population of P(t) individuals, for generation t. Each of these individuals contain a potential solution to the posed problem that need to be evaluated and its fitness measured. Some of these individuals are randomly selected to undergo a stochastic transformation and become new individuals (genetic operation). Likewise natural genetics, this transformation can be a mutation, which creates new individuals by making changes in a single individual, or crossover, which creates new individuals by combining parts from two others. The offspring from this process, the new individuals C(t), are evaluated and its fitness measured. A new population is created after selecting the more fit individuals from the parent and the offspring population. In the end, after several generations, the algorithm converges to the best individual, which is a possible optimal or sub-optimal solution to the problem, as stated by Gen and Cheng [13].

The genetic algorithm components need to be carefully addressed in order to provide a good search space and exploit the best solution. A good balance between exploration and exploitation is a must for complex and real-world problems.

In a mesh-free discretization, the size of the compact support, where nodal shape functions are defined, and the size of the domain of integration, where the nodal stiffness matrix of the numerical model is computed, must be conveniently defined in any application, since their values strongly affect the performance of the numerical solution. Therefore, the values of the size of the compact support and the values of the size of the local integration domain, are optimized in this paper. They are IV. Multi-Objective Optimization of the Mesh-free Model defined, respectively in Equations ( 10) and (11) which show that the accuracy of a mesh-free numerical application can be controlled through a proper specification of the discretization parameters ? s and ? q . Therefore, parameters ? s and ? q are both set as design variables of the multi-objective optimization process, in order to be automatically defined with optimal values. Additionally, in order to facilitate and automate the pre-processing phase of the mesh-free modeling, the nodal distribution need to be addressed. Therefore, for a bi-dimensional problem, the number of divisions in x and y direction are chosen as design variables. When the number of divisions in both directions are provided, the mesh-free numerical model, can define the nodal coordinates and distribute the nodes along the problem domain and boundary, including crack nodes. For this case, only regular nodal distributions are considered.

The use of efficient objective functions condition the overall performance of the multi-objective optimization process. The objective functions force, through meshfree numerical simulation, the minimum total mechanical energy of the structure and the conditioning of the final system of algebraic equations which, consequently enforce the solution accuracy of the mesh-free model. Note that solution errors of the mesh-free model are not included in any of the objective functions which, therefore are quite general and do not depend on any analytical solution.

The standing challenge in the application of numerical simulations in the optimization process is the accurate evaluation of objective functions which obviously is dependent upon the automatically generated mesh-free discretization. Since multiple iterations are required during the optimization process, it is necessary to maintain a balance between efficiency and accuracy through constraints of the design variables.

The definition of this objective function results from the features of the parameter ? s in combination with ? q . Considering a body with the actual elastic field in any state, the strain energy U , and the potential energy P , of external forces, respectively given by
U = ? 1 2 ? T ? d?(42)
and
P = ? Î?"t t T u dÎ?",(43)
can be used to obtain the total potential energy T . The work theorem, when applied to the global domain of the body, for the actual elastic field settled in the body, leads to P = ?2U and therefore T = ?U , as well as T = P/2. These results show that the minimum value of the total potential energy of the body corresponds to a minimum value of the potential energy P or a maximum value of the strain energy U . The energy can be measure both by strain energy U or potential energy P , although evaluation of U is computationally more expensive, since requires the computation of the stress field for all nodal values and derivatives of shape functions. Therefore, the potential energy P is used instead, since it requires the evaluation of displacement fields only at static boundary nodes, the ones with no-null ap- plied loads, and does not require the computation of derivatives of shape functions, which is computationally efficient in comparison.

Hence, the objective function can be defined with the structural compliance C,
as C = 1 2 Î?"t t T u dÎ?" = ? 1 2 P.(44)
Consequently, the minimum value of the potential energy P corresponds to a maximum value of ?C that is equivalent to a minimum value of C.

The numerical problem optimization aims to minimize the objective function using the mesh-free numerical model, by finding optimum values for the design variables, in this case the geometrical parameters ? s , ? q and the nodal distribution, also satisfying the problem constraints.

The mathematical formulation of the multi-objective optimization scheme for linear elastic fracture mechanics problems is as follows minimize C(? s , ? q , n x , n y )

CPU time(? s , ? q , n x , n y )
subject to e(? s ) = ? s min ? ? s ? ? s max e(? q ) = ? q min ? ? q ? ? q max e(n x ) = n x min ? n x ? n x max e(n y ) = n y min ? n y ? n y max
where ? s = (? s1 , ? s2 , ..., ? sn ) ? ? s ? q = (? q 1 , ? q 2 , ..., ? q n ) ? ? q n x = (n x1 , n x2 , ..., n xn ) ? (x)
n y = (n y 1 , n y 2 , ..., n y n ) ? (y) ,(45)
in which C is the structural compliance, CPU time is the time required to generate and solve the global system of algebraic equations; ? s min /? q min and ? s max /? q max denote the minimum and the maximum allowable limits for the mesh free discretization parameters ? s and ? q , respectively. n min /n max denote the minimum and the maximum geometrical values for the number of divisions on both directions (x and y), limited by the geometrical constraint of the problem, for a regular nodal discretization of the posed problem. Therefore, the variable n also determine the total number of nodes for the problem and node coordinates, automatically defined for a regular nodal distribution. On this multi-objective optimization, the fitness function, that is the routine containing the mesh-free algorithm, define scalar values for ? s , ? q , n x and n y , yielding different objective function outputs. Since there are two objective functions, the Pareto front will be the final result of the optimization, which will provide non-dominant solutions.

The ILMF is the only mesh free method implemented in this paper, but this process can be easily applied to any desired local mesh free method. The whole optimization process is summarized in the flowchart presented in Figure 3. This section presents numerical results to demonstrate the accuracy and efficiency of the mesh-free numerical method with optimization, through different linear fracture mechanics problems previously presented by Oliveira and Portela [26].


# d) Algorithm Formulation


# V. Numerical Results

For a regular mesh-free discretization of n x x n y nodes, the size of the local support ? s and the size of the local integration domain ? q , are respectively parameters ? s and ? q . Good results can be obtained with a mesh-free model if r ?s , r ?q and the arrangement of nodes are properly refined.

Usually, these parameters and the nodal distribution are heuristically defined. One key advantage of the ILMF modeling process is that it can provide appropriate Flowchart of the multi-objective optimization scheme for mesh-free numerical methods. values for ? s and ? q using genetic algorithms, as initially presented by Santana et al. [31], which greatly improves the model accuracy. Additionally, this work also optimize the nodal distribution, resulting in a fully automated optimization routine for the entire pre-processing phase of traditional numerical methods, which is the definition of the mesh.

Three cases of edge-cracked square plates, respectively under mode-I, mode-II and mixed-mode deformation are considered.

The discontinuity generated by the presence of the crack requires a special treatment in order to be carried out in this non-convex domain. Therefore, the crack faces are modeled with two lines of overlapping nodes, where the MLS approximation is acting only in their respective influence size, while the crack tip is modeled with one node that can influence both sides of the crack. The visibility criterion is implemented around the crack during the definition of the compact support of each node. Hence, the compact support and the local integration domain of each node of the crack faces are defined as in the case of a traction-free boundary node.

The size of the local integration domain of the crack tip node is defined as ? k = ? q /2, to ensure the local aspect of the discretization of the crack. The computation of matrices g and f , of the Williams' singular solution at each crack tip is carried out with Gaussian quadrature, with a single integration point.

The results obtained with the ILMF using the multi-objective optimization are compared with the results originally obtained by Oliveira and Portela [26], without optimization, and by Portela and Aliabadi [28], using the DBEM with the J-integral (J-DBEM) technique, which proved to be a very accurate method. The DBEM modeling strategy considers piecewise-straight cracks which are discretized with straight discontinuous quadratic boundary elements. Continuous quadratic boundary elements are used along the remaining boundaries of the problem, except at the intersection between a crack and an edge, where semi-discontinuous boundary elements are used on the edge. Self-point discontinuous boundary elements are integrated analytically, while Gaussian quadrature, with sub-element integration, is carried out for the remaining integrations.

The GA is set to minimize the Compliance C and CPU time or computational effort, chosen as objective functions for this optimization process. The design variables of the optimization, the number of nodes and the node coordinate are defined within the problem geometry, and the parameters ? s = 1.5 ? 10 and ? q = 0.1 ? 0.9, are defined as continuous in the intervals. Only the major computational cost that is the cost of generating and solving the global system of algebraic equations, was measured.

On this optimization scheme, the initial population is randomly generated according to the predefined population size of 25 individuals. Then, the fitness function is calculated for each member of the population and scaled using a rank process, which is used later in the selection process. The reproduction operator is implemented based on a tournament selection. Both mutation and crossover are constraint dependent. The genetic algorithm described above generates a stochastic values sequence of design variables which are evaluated through the objective functions. Finally, the optimization process is terminated if the number of generations exceeds the predefined maximum number, which is selected as 150 in this scheme, or if the average change in fitness function is less than 1 × 10 ?6 . The improved accuracy of the optimization process can be clearly seen on this benchmark problem, regardless of the loading.

A square edge-cracked plate, represented in Figure 4, is considered for the first analysis. The plate, with crack length a, width w and height h = w/2, is loaded
a W 2 / W 2 / W
Square plate with a single edge crack under mode-I loading (h/w = 0.5). by a uniform traction t = ?, applied symmetrically at the ends. All the results presented are for h/w = 0.5, to be compared with the highly accurate values introduced by Civelek and Erdogan [7]. Therefore, five cases were considered, with a/w = 0.2, 0.3, 0.4, 0.5 and 0.6.

The ILMF model was applied with rectangular local domains of integration, with discretization parameters and nodal configuration automatically defined though GA optimization. The MLS approximation considered a first-order polynomial basis with quartic spline weighting function. It is important to highlight that all nodal distributions were performed without considering any refinement of the discretization around the crack tip, always with regular distributions.

Figure 5 show the Pareto front obtained from the optimization process, containing all feasible solutions for the posed problem. From the frontier solutions, The multi-objective Pareto front of the square plate with a single edge crack under mode-I loading, for a/w = 0.5; ILMF with the automatic parameters optimization routine. a set of solutions were selected and the results presented in Figure 6 and Table 1; where it can be seen that the optimization lead to accurate results for all points The multi-objective Pareto front results of selected feasible solutions for the square plate with a single edge crack under mode-I loading, for a/w = 0.5. in the Pareto front, with minimum values for compliance and SIF close to reference values. For this case, ? s greatly varies depending on the nodal distribution, but the best values for ? q are usually closer to 0.5. Table 2 show the results obtained for different a/w, where ILMF represents the values obtained in Oliveira and Portela [26], ILMF + represents the values presented in this work using the optimization routine, Portela and Aliabadi [28] represents the values obtained with Table 1: The multi-objective Pareto front of selected feasible solutions for a square plate with a single edge crack under mode-I loading (a/w = 0.5).


# i. Mode-I Loading
Index
Square plate with a single edge crack under mode-I loading.
K I /(t ? ?a)
% Error a/w ILMF ILMF + J-DBEM [28] Reference [7] ILMF + J-DBEM 0.2 1.520  [7]. In this analysis, the SIF values of the mode-II are always below 10 ?7 , since this is a mode-I loading crack problem.

The results highlight the accuracy of ILMF, which was further improved after the optimization process, always very close to reference values and J-DBEM. Even for similar nodal distributions as originally conceived, like index 1 of the Pareto front of selected solutions, the optimization of ? s was enough to improve the overall accuracy. The nodal distribution obtained by the GA optimization scheme and the respective deformed configuration of the plate is schematically represented in Figures 7 and 8. Regular nodal distribution resulting from the optimization scheme, with a regular nodal distribution of 20 × 18 = 360 nodes and additional overlapping nodes on the crack faces, for a/w = 0.5, under mode-II loading. The red line represents the crack faces.


# Table 2:

A square edge-cracked plate, with ratio between the height and the width of the plate as h/w = 0.5, schematically represented in Figure 9, is considered for this analysis.

The plate is loaded with a uniform traction t, parallel to the crack of length a and is applied anti-symmetrically on the sides which corresponds to a mode-II loading. There are no published benchmark results due to the complexity of the problem and, therefore, they are compared with the results obtained with the J-DBEM, using the software [28].  For this problem, five cases were considered, with corresponding ratios of a/w = 0.2, 0.3, 0.4, 0.5 and 0.6. Rectangular local domains of integration, firstorder polynomial basis and quartic spline weighting function are considered on the ILMF model. Like the previous problem, the regular nodal configuration and discretization parameters are automatically defined though GA optimization, without any special refinement around the crack tip.

The results obtained for this multi-objective optimization process are presented in Figure 10, Figure 11 and Table 3; with all point in the Pareto front leading to The multi-objective Pareto front of the square plate with a single edge crack under mode-II loading, for a/w = 0.5; ILMF with the automatic parameters optimization routine.

The multi-objective Pareto front results of selected feasible solutions for the square plate with a single edge crack under mode-II loading (a/w = 0.5).   accurate results and fast computations. It can be seen that SIF values are close to each other, depending on the minimum value of the compliance indicator. For this optimization, both ? s and ? q values are close to each other due to the similarity between nodal distributions obtained. Once more, ? q ? 0.5 and, for this case, ? s = 1.4 ? 4. The multi-objective Pareto front of selected feasible solutions for a square plate with a single edge crack under mode-II loading (a/w = 0.5).


# Index CPU Time

For different values of a/w, Table 4 show the results obtained from the anal-Square plate with a single edge crack under mode-II loading.  [28]. In this problem, the SIF values obtained for the mode-I are always below 10 ?3 , since this is a mode-II crack problem.
K I /(t?
Even though this problem is highly complex, accurate values were obtained after the optimization process, improving the previous results. The nodal distribution obtained by the GA optimization scheme and the respective deformed configuration of the plate is schematically represented in Figures 12 and 13.

Consider now a plate with an edge slant crack, as represented in Figure 14 schematically, in mixed-mode deformation. The length of the crack is denoted by a, the width and height of the plate is denoted by w. The plate is loaded by a uniform traction t = ?, applied symmetrically at the ends.

iii. Mixed-Mode Loading Table 4:   tios of a/w = 0.2, 0.4 and 0.6, for ? = 30 ? and two cases, with corresponding ratios of a/w = 0.2 and 0.4, for ? = 60 ? , as originally presented by Murakami [24]. For MLS approximation of the elastic field, a first-order polynomial basis and a quartic spline weighting function were considered, along with rectangular local domains to perform the numerical integration of the ILMF model. The regular nodal configuration and discretization parameters are automatically defined though GA optimization, without any special refinement around the crack tip.

All the results presented are compared with the accurate values provided by Murakami [24] and Portela and Aliabadi [28]. The multi-objective optimization process resulted in the Pareto front of Figure 15, where all feasible solutions are presented. The selected multi-objective optimization are presented in Figure 16 and Tables 5 and 6; where a good accuracy can be seen related to minimum values The multi-objective Pareto front results of selected feasible solutions for the square plate with a single edge crack under mixed-mode loading, for a/w = 0.4 and ? = 30   The multi-objective Pareto front of selected feasible solutions for the square plate with an edge slant crack under mixed-mode loading (a/w = 0.4).

The multi-objective Pareto front results of selected feasible solutions for the square plate with a single edge crack under mixed-mode loading, for a/w = 0.4 and ? = 60 from the values obtained with the J-integral implemented in the DBEM, provided by Portela and Aliabadi [28], and Murakami [24]. The nodal distribution obtained by the GA optimization scheme and the respective deformed configuration of the plate is schematically represented in Figures 17 and 18.

The compliance proved to be an efficient objective function for linear elastic fracture mechanics problems and complement the already efficient SST implementation, without any refinement around the crack tip due to the regularized stress field.

Square plate with a single edge crack under mixed-mode loading, for ? = 60 ? and K I /(t ? ?a).

K I /(t ? ?a) % Error a/w ILMF ILMF + J-DBEM [28] Reference [24]     Deformed configuration of the plate resulting from the optimization scheme, for a/w = 0.6, under mixed-mode loading.

The ILMF local mesh free numerical method, implemented with SST, was improved through an optimization scheme that automatically define the nodal distribution and the discretization parameters, for solving two-dimensional problems of the linear elastic fracture mechanics.

The MLS and reduced numerical integrations are considered in the discretization of the elastic field, using a node-by-node process to generate the global system of equilibrium equations, which is very efficient and prone to parallel processing. Also, the reduced integration reduce the stiffness associated with local nodes, leading to an increase in the overall accuracy, without the well-known instabilities associated with the process.

The SST implemented for linear elastic fracture mechanics applications performs a regularization of the stress field, introducing the SIF as additional primary unknowns of the problem. As a consequence, the analysis does not require refined nodal distributions around crack tips, in contrast to other numerical methods. The numerical results are evidence of the efficiency of the modeling strategy, since accurate results were obtained for edge-cracked square plates under mode-I, mode-II and mixed-mode, always without any refinement around the crack tip and relatively small nodal distributions, automatically obtained by the optimization algorithm.

Historically, the nodal distribution, the size of the compact support and the size of the local integration domain are heuristically defined and need to be addressed for every mesh-free application. The ILMF model has the capability of the automatic definition of the discretization parameters and the nodal distribution, through a multi-objective optimization process, based on GA.

The definition of the objective function as a profound impact on the behavior of the optimization process and need to be carefully defined. In this paper, an appropriate objective function is derived from the classical structural theorem of the minimum total potential energy, carried out only at static boundary nodes that does not require the computation of derivatives of shape functions. Therefore, the optimization scheme is computationally very efficient and as the additional benefit of not requiring the analytical solution to be performed.


# VI. Conclusions

The results obtained with the optimization algorithm are in agreement with those of the reference values, where low compliance values are associated with accurate SIF values, as expected. This result show that the local Pareto-optimal is always quite close to the global Pareto-optimal solutions, which is always desirable from a computational point of view. The structural compliance objective function effectively optimized the discretization parameters and the nodal distribution, properly defining these geometrical properties with fast computations and without any user input.

This paper show that mesh-free methods, along with optimization processes, could provide stable and accurate solutions for fracture mechanics problems with minimal user input, contributing to a mainstream use of mesh-free numerical methods in the near future.
























3(s) Compliance K I /(t ? ?a)? s? qNodes10.3949.87E-040.2642.662 0.51123720.2410.2040.3151.485 0.50129930.967-0.1150.2813.903 0.50325140.8860.7910.2313.885 0.828359
5Ind. CPU T.(s)CK I /(t? ?a) K II /(t ? ?a)? s? qN.13.082.48E-041.6670.5055.821 0.617 30126.614-4.68E-021.5180.4559.769 0.612 45530.3294.46E-041.6450.4942.753 0.916 20345.744-1.34E-041.9370.5876.56 0.499 33151.4372.54E-041.5790.4794.25 0.501 267for compliance.? .
67Ind. CPU T. (s)CK I /(t? ?a) K II /(t? ?a)? s? qN.10.9094.61E-040.6490.4544.761 0.503 20321.9156.84E-040.7010.4918.58 0.501 16538.2172.29E-040.7130.4997.955 0.498 19940.4017.08E-040.5110.4823.659 0.495 199compared to other examples or ? Square plate with a single edge crack under mixed-mode loading, for ? ?a). ? = 30 ? and K I /(t K I /(t ? ?a) % Errora/w ILMF ILMF + J-DBEM [28] Reference [24] ILMF + J-DBEM0.21.1641.11.0821.1000.0090.0160.41.513 1.5791.5451.5500.0180.0030.62.732 2.7432.5722.5500.080.009? .
ILMF + J-DBEM0.20.543 0.5200.4950.5000.0390.0100.40.603 0.6040.5920.6000.0050.013K II /(t? ?a)% Errora/w ILMF ILMF + J-DBEM [28] Reference [24] ILMF + J-DBEM0.20.327 0.3730.3560.3600.0350.0110.40.439 0.4220.4130.4200.0050.01700.250.50.751
8910
			© 2021 Global Journals
		
		
			

				


* 
	
		The meshless local petrov-galerkin (mlpg) method: A simple and less-costly alternative to the finite element and boundary element methods
		
			SNAtluri
		
		
			SShen
		
	
	
		CMES: Computer Modeling in Engineering and Sciences
		
			3
			1
			
			2002
		
	




* 
	
		A new meshless local petrov-galerkin (mlpg) approach in computational mechanics
		
			SNAtluri
		
		
			TZhu
		
	
	
		Computational Mechanics
		
			22
			2
			
			1998
		
	




* 
	
		Optimization of meshless local petrov-galerkin using genetic algorithm for 3d elastostatic problems
		
			ABagheri
		
		
			REhsany
		
		
			MJMahmoodabadi
		
		
			GHBaradaran
		
	
	
		International Journal of Engineering
		
			24
			2
			
			2011
		
	




* 
	
		Optimal pareto parametric analysis of two dimensional steady-state heat conduction problems by mlpg method
		
			GHBaradaran
		
		
			MJMahmoodabadi
		
	
	
		International Journal of Engineering
		
			22
			4
			
			2009
		
	




* 
	
		Cracked plate analysis with the dual boundary element method and williams' eigenexpansion
		
			JCaicedo
		
		
			APortela
		
	
	
		J. of Engineering Analysis with Boundary Elements
		
			52
			
			2015
		
	




* 
	
		Meshfree methods: Progress made after 20 years
		
			JChen
		
		
			MHillman
		
		
			SChi
		
	
	
		Journal of Engineering Mechanics
		
			143
			4
			4017001
			2017
		
	




* 
	
		Crack problems for a rectangular plate and an infinite strip
		
			MBCivelek
		
		
			FErdogan
		
	
	
		International Journal of Fracture
		
			19
			
			1982
		
	




* 
	
		Multiphysics topology optimization of heat transfer and fluid flow systems
		
			EDede
		
	
	
		Proceedings of the COMSOL Conference
				the COMSOL ConferenceBoston
		
			COM-SOL
			2009
		
	




* 
	
		Isogeometric analysis for topology optimization with a phase field model
		
			LDede
		
		
			MJBorden
		
		
			TJ RHughes
		
	
	
		Archives of Computational Methods in Engineering
		
			19
			
			2012
		
	




* 
	
		References Références Referencias
		
	




* 
	
		Nodal cosine sine material interpolation in multi objective topology optimization with the global criteria method for linear elasto static, heat transfer, potential flow and binary cross entropy sharpening
		
			MDenk
		
		
			KRother
		
		
			MZinßer
		
		
			CPetroll
		
		
			KPaetzold
		
	
	
		Proceedings of the 23th International Conference on Engineering Design
				the 23th International Conference on Engineering Design
		
			Design Society
			2021
			
		
	




* 
	
		Nodal cosine sine material interpolation in multi objective topology optimization with the global criteria method for linear elasto static
		
			MDenk
		
		
			KRother
		
		
			MZinßer
		
		
			CPetroll
		
		
			KPaetzold
		
		
			2021
		
	
	heat transfer, potential flow and binary




* 
	
		Adaptive refinement in the meshless finite volume method for elasticity problems
		
			MEbrahimnejad
		
		
			NFallah
		
		
			ARKhoei
		
	
	
		Computers & Mathematics with Applications
		
			69
			12
			
			2015
		
	




* 
	
		Genetic Algorithms and Engineering Optimization
		
			MGen
		
		
			RCheng
		
		
			2000
			Wiley
		
	




* 
	
		Topology optimization of heat conduction problems using the finite volume method. Structural and Multidisciplinary Optimization
		
			AGersborg-Hansen
		
		
			MPBendsoe
		
		
			OSigmund
		
		
			2006
			31
			
		
	




* 
	
		Topology optimization of front metallization patterns for solar cells. Structural and Multidisciplinary Optimization
		
			DKGupta
		
		
			MLangelaar
		
		
			MBarink
		
		
			FVan Keulen
		
		
			2015
			51
			
		
	




* 
	
		Adaptation in Natural and Artificial Systems
		
			JHHolland
		
		
			1975
			MIT press
		
	




* 
	
		
			AHuerta
		
		
			TBelytschko
		
		
			SFernández-Méndez
		
		
			TRabczuk
		
		
			XZhuang
		
		
			MArroyo
		
		Meshfree methods. Encyclopedia of Computational Mechanics
				
			ErwinStein
			RenéDe Borst
			ThomasJ RHughes
		
		
			Wiley
			2017
		
	
	Second Edition




* 
	
		Multiple Objectives Decision Making-Methods and Applications
		
			CLHwang
		
		
			AS MMasud
		
		
			1979
			Springer
		
	




* 
	
		Application of genetic algorithms to lubrication pump stacking design
		
			VKelner
		
		
			OLeonard
		
	
	
		Journal of Computational and Applied Mathematics
		
			168
			1
			
			2004
		
	




* 
	
		Multiobjective evolutionary structural optimization using combined static/dynamic control parameters
		
			WYKim
		
		
			RVGrandhi
		
		
			MHaney
		
	
	
		AIAA J
		
			44
			
			2006
		
	




* 
	
		A local point interpolation method for stress analysis of two-dimensional solids
		
			GRLiu
		
		
			YTGu
		
	
	
		Structural Engineering and Mechanics
		
			11
			2
			
			2001
		
	




* 
	
		Point interpolation method based on local residual formulation using radial basis functions
		
			GRLiu
		
		
			LYan
		
		
			JGWang
		
		
			YTGu
		
	
	
		Structural Engineering and Mechanics
		
			14
			
			2002
		
	




* 
	
		Genetic algorithms for modelling and optimisation
		
			JMccall
		
	
	
		Journal of Computational and Applied Mathematics
		
			184
			1
			
			2005
		
	




* 
	
		
			YMurakami
		
		Stress intensity factors handbook
				
			Pergamon Press
			1987
			2
		
	
	1st edition




* 
	
		Weak form collocation -a local meshless method in linear elasticity
		
			TOliveira
		
		
			APortela
		
	
	
		Engineering Analysis with Boundary Elements
		
			73
			1
			
			2016
		
	




* 
	
		A local mesh free method with the singularity subtraction technique
		
			TOliveira
		
		
			APortela
		
	
	
		Engineering Analysis with Boundary Elements
		
			104
			
			2019
		
	




* 
	
		A local mesh free method for linear elasticity and fracture mechanics. Engineering Analysis with Boundary Elements
		
			TOliveira
		
		
			WVélez
		
		
			ESantana
		
		
			TAraújo
		
		
			FMendonça
		
		
			APortela
		
		
			2019
			101
			
		
	




* 
	
		Crack Growth Analysis Using Boundary Elements -Software
		
			APortela
		
		
			MHAliabadi
		
	
	
		Proceedings of the 23th International Conference on Engineering Design
				the 23th International Conference on Engineering DesignSouthampton, UK and Boston, USA
		
			Design Society
			1993
			
		
	




* 
	
		Multiobjective Optimization: Behavioral and Computational Considerations
		
			JLRinguest
		
		
			1992
			Kluwer
		
	




* 
	
		A local mesh free numerical method with automatic parameter optimization
		
			ESantana
		
		
			TOliveira
		
		
			WVélez
		
		
			TAraújo
		
		
			FMartins
		
		
			APortela
		
	
	
		Engineering Analysis with Boundary Elements
		
			113
			
			2020
		
	




* 
	
		Theory of Multiobjective Optimization
		
			YSawaragi
		
		
			HNakayama
		
		
			TTanino
		
		
			1985
			Academic Press
		
	




* 
	
		Multiple Criteria Optimization: Theory, Computation, and Application
		
			RESteuer
		
		
			1986
			Wiley
		
	




* 
	
		Stress singularities resulting from various boundary conditions in angular corners of plates in extension
		
			MLWilliams
		
	
	
		Journal of Applied Mechanics
		
			
			1952
		
	




* 
	
		No free lunch theorems for optimization
		
			DHWolpert
		
		
			WGMacready
		
	
	
		Transactions on Evolutionary Computation
		
			1
			1
			
			1997
		
	




* 
	
		A comparative evaluation of genetic and gradient-based algorithms applied to aerodynamic optimization
		
			DavidWZingg
		
		
			MarianNemec
		
		
			ThomasHPulliam
		
	
	
		European Journal of Computational Mechanics
		
			17
			1-2
			
			2008
		
	




* 
	
		Topology-optimized 4d printing of a soft actuator
		
			AZolfagharian
		
		
			MDenk
		
		
			MBodaghi
		
		
			AZKouzani
		
		
			AKaynak
		
	
	
		Acta Mechanica Solida Sinica
		
			33
			
			2020
		
	




* 
	
		Multicriterion evolutionary structural optimization using the weighting and the global criterion methods
		
			KProos
		
		
			GSteven
		
		
			OQuerin
		
		
			YXie
		
	
	
		AIAA J
		
			39
			2006-2012, 2001