Comparative Study of Peridynamics and Finite Element Method for Practical Modeling of Cracks in Topology Optimization

: Recently, topology optimization of structures with cracks becomes an important topic for avoiding manufacturing defects at the design stage. This paper presents a comprehensive comparative study of peridynamics-based topology optimization method (PD-TO) and classical ﬁnite element topology optimization approach (FEM-TO) for designing lightweight structures with/without cracks. Peridynamics (PD) is a robust and accurate non-local theory that can overcome various difﬁculties of classical continuum mechanics for dealing with crack modeling and its propagation analysis. To implement the PD-TO in this study, bond-based approach is coupled with optimality criteria method. This methodology is applicable to topology optimization of structures with any symmetric/asymmetric distribution of cracks under general boundary conditions. For comparison, optimality criteria approach is also employed in the FEM-TO process, and then topology optimization of four different structures with/without cracks are investigated. After that, strain energy and displacement results are compared between PD-TO and FEM-TO methods. For design domain without cracks, it is observed that PD and FEM algorithms provide very close optimum topologies with a negligibly small percent difference in the results. After this validation step, each case study is solved by integrating the cracks in the design domain as well. According to the simulation results, PD-TO always provides a lower strain energy than FEM-TO for optimum topology of cracked structures. In addition, the PD-TO methodology ensures a better design of stiffer supports in the areas of cracks as compared to FEM-TO. Furthermore, in the ﬁnal case study, an intended crack with a symmetrically designed size and location is embedded in the design domain to minimize the strain energy of optimum topology through PD-TO analysis. It is demonstrated that hot-spot strain/stress regions of the pristine structure are the most effective areas to locate the designed cracks for effective redistribution of strain/stress during topology optimization.


Introduction
Many structural systems contain multiple elements that need to be designed as light as possible for various engineering applications such as aerospace [1,2], automotive [3,4], and other structural design industries [5][6][7]. Prior to manufacturing these lightweight structures, their designs are commonly performed using density, shape and topology optimization(TO) [8][9][10] methods. Topology optimization finds out the most effective and optimal material deposition/distribution on a predefined design domain by optimizing a given objective function subjected to relevant constraints.
Earlier TO methods involved ill-posed mathematical algorithms that may lead to in-stable results or occurrence of checkerboard pattern in final topology. Bendsøe et al. [8] solved this problem by proposing a homogenized TO method. Later, more popular methods including the Solid Isotropic Material with Penalization (SIMP) [11,12] technique and the Evolutionary Structural Optimization (ESO) technique [13] were developed to serve highly efficient and accurate TO strategies. For instance, the SIMP method was used in [14] for designing multi-material compliant structures. In addition, the ESO method was employed for multi-material interpolation [15], whereby a combination of constraints and a heuristic design variables updating scheme were considered to perform the compliance minimization. Afterwards, bidirectional ESO (BESO) method was implemented by Querin et al. [16] to check all possible directions, including material removal and also addition of material to elements according to their high-stress levels.
In general, TO methods can be classified based on discrete or continuous design variables. For instance, discrete-based approaches were implemented using genetic algorithms [17,18], simulated annealing optimization [19], and non-gradient topology optimization (NGTO) [20] to solve TO problems. In addition, continuous density based methods are employed together with the SIMP algorithm. Essentially, the continuous methods can be categorized as two main approaches. Its first kind is the proportional approach (PO) [21], in which the density of the elements is assigned according the value of objective function in the previous iteration. Apart from PO approach, two decades ago, Sigmund [22] proposed the optimality criteria approach (OC) for effectively solving TO problems using continuous density variables. This approach attempts to satisfy a set of analytically obtained criteria instead of directly optimizing the objective.
During TO, a proper numerical method is needed to perform accurate structural analysis. To this end, one of the most common technique is classical continuum mechanics (CCM) formulations that use finite element method (FEM) [23,24]. Since CCM formulation is established according to surface traction of locally interacting material points/particles, it mostly deals with a structure as a whole [25,26], and therefore the domain remains continuous as the structure deforms. For deriving mathematical relations according to mass, momentum and energy principles in CCM, particles interactions are considered between a particle and its nearest neighbor. However, such assumptions pose a modeling/analysis limitation to the CCM applications for structures including damage, discontinuity, internal feature or defect. In other words, CCM is incapable of dealing with discontinuity since it use spatial partial differential equations. In the last decades, various research efforts have been dedicated to overcome the limitations of CCM for handling cracked structures such as Linear Elastic Fracture Mechanics (LFEM) [27,28] and eXtended Finite Element Method (XFEM) [29,30]. LFEM dealt with crack near its tip as a singular point and dealt with critical energy release rate (artifact manipulation of equation). However, for a complex problem with multiple interacting cracks, this method still becomes complex since it causes the requirement for redefining the body at each step. Thus, solving these structures with traditional finite element is very complicated. Furthermore, XFEM overcomes some of these issues as it uses partition of unity in domain regarding property of finite elements [31]. Later, the near crack tip elements are partially enriched [32] in the XFEM. Despite these assumptions, XFEM solutions somewhat become mostly inaccurate especially in the blending regions.
Another approach for accurately performing structural analysis during TO is non-local continuum theories. These approaches can simply eliminate inefficiency and computational complexity of CCM formulations for modeling dynamic/static problems of cracked structures. One of most common non-local theories that have been extensively studied in the last decade is the reformulation of continuum mechanics, referred to as Peridynamics (PD) [33,34]. Unlike CCM, the PD uses integrodifferential equations so it doesn't contain any spatial derivatives. In this regard, the material does not necessarily require to remain continuous during the simulations. As such, PD becomes ideal approach to deal with discontinuities such as cracks and to predict crack initiation and also crack propagation without adding complexity to the system of equation. Overall, these features make PD a viable tool to nucleate damage in the domain of structure without adding any singularity.
Since its first publication, Silling [35] firstly used PD to simulate complex crack growth in plate structures. Then, the PD have been widely applied to fracture mechanics simulation of various material and structural systems [36][37][38][39][40]. Further studies such as [41,42] used PD to solve multi-scale structures failure and damage prediction of different materials in composites and layered heterogeneous materials [43]. Recently, Basoglu [44] introduced the PD for integrating micro-crack in structure. PD is also used for predicting different type of fracture modes [45] on various structures including concrete [46], anisotropic [47] and functionally graded materials [48]. There are very few studies included direct utilization of PD for TO of cracked structures. The first study in literature was performed by Kefal et al. [49], who integrated bond-based PD with BESO optimization scheme to perform TO of structures with discontinuities. Bond-based PD can be considered as a non-local method accommodating symmetric interactions between particles. Later, peridynamics based topology optimization (PD-TO) approach was extended to [50] gradient-based optimization algorithms for complex engineering problems involving cracks. Most recently, the PD-TO algorithm was also utilized for multi material topology optimization [51].
To the best of authors' knowledge, there is no research study on the extensive comparison of PD-TO capabilities for topology optimization of cracked structures with respect to FEM-TO. In fact, the available studies ( [49,50]) only make consecutive comparison between FEM and PD for performing topology optimization of structures without cracks. However, there is a need for justifying why TO algorithms should substitute peridynamics for conventional FEM approach in topology optimization of cracked structures. This important question can only be addressed if the capabilities of the existing FEM-based optimization algorithms and PD-TO are utilized and compared on the same design domain having an initial crack. The main and novel aim of the present study is to resolve this challenging problem by performing comprehensive simulations on scrupulously selected TO case studies with single/multiple cracks. Here, to be able to model the cracks practically in the FEM analysis, the finite elements encountered in the cracked regions are removed from the design domain during the optimization. On the other hand, the PD bonds are broken on these same regions of the crack line during the PD-TO simulations. Apart from this aspect, the present study has another crucial novelty in terms of proposing a useful computational strategy for designing lighter-weight structures via integration of possible cracks at the design stage. In other words, this study proposes the utilization of PD-TO for locating cracks in such way that the final topology of cracked structure becomes better (stiffer) than the pristine one (i.e., having no cracks). Such an approach can be suitable for modeling, analysis, and optimization steps using PD since the other methods may be incapable of modeling crack tip singularities during TO steps. Thus, this study also sheds a light to enhancement of the topology optimization via PD and cracks and increases the technological readiness level of the PD-TO methodology by proposing a viable computational strategy for the first time in literature. The advancements studied herein will be beneficial to enable the peridynamics and topology optimization community for designing superior structures with higher stiffness-to-weight ratio.
This paper organized as follows. First, peridynamic theory is briefly explained in Section 2.1. Then, mathematical formulation of the optimality criteria approach is presented together with the description of filtering concept in Section 2.2. Afterwards, four case studies are investigated in Section 3. Except the third benchmark problem (a bracket which is studied for effect of geometrical complexity), all case studies are considered to have symmetrical geometry and boundary condition to reflect the effect of symmetry in the topology optimization of structures with/without cracks. For each case, PD-TO results of the structures without cracks are compared to those predicted by FEM-TO. After that, different cracks are included in the design domains of the structures to compare TO capability of FEM versus PD for the cracked structures. For comparison between local and non-local approaches, the optimized topologies are assessed by evaluating the strain energy and displacement results. For the last case study, effect of crack in topology optimization is also studied to decrease strain energy by integrating a crack in a stress-concentrated region. Finally, concluding remarks of the study are provided in Section 4.

Peridynamic
In this section, we will summarize the mathematical formulation of Peridynamicsbased topology optimization. PD theory is a non-local continuum theory, where the structural media is represented by the ensemble of infinitely small volumes, so-called material points (particles). Unlike local continuum mechanics, each particle interact with all particles within its support region, named as Horizon, in PD. Figure 1 shows a set of particles in a region β and horizon of particle i with size of δ. For each material point i located at position x , the horizon can be described as: According to the bond-based theory, the interaction of particles are defined in symmetrical pairwise manner, thus the bond forces acting on particle x and x have equal magnitude but with opposite sign f and f , f = − f .
The pairwise body force density f between x and x is function of relative position vector ξ = x − x and displacement vector η = u − u as f(u − u, x − x, t). According to [34], the equation of motion for sets of particle is defined as : where ρ(x), b(x, t) andü(x, t) are density, body force and acceleration of particle x, respectfully. Relative position vector of two material points in deformed configuration can be written as y − y = (u − u) − (x − x). Accordingly the stretch between these two material points can be defined as: Utilizing the stretch and peridynamic material parameters, the force density vector can be expressed along the bond direction in the deformed configuration as: where c is the bond constant that can be calculated through equating the strain energy densities obtained from PD and CCM theories for an isotropic material. For instance, this material parameter can be explicitly given under plane stress condition for 2-D domain as [33]: where E is young modulus of isotropic material and h is thickness of the domain. Moreover, micro-potential energy of each bond can be calculated using linear force and displacement relation of the bond as: Integrating this micro-potential for each bond in the horizon of particle x, the strain energy density of each particle can be obtained as: To be able to obtain total strain energy of domain β, one can integrate the strain energy density over the full field domain as: Here U is the total strain energy that is used as a compliance in the topology optimization process.
To include local cracks into peridynamic simulation, ratio of the number of eliminated interactions to the total number of initial interactions is defined as a weighted function φ: where µ is damage initiation value which is a history-dependent scalar-valued function that can be defined as : where s c is the critical stretch value. This value is dependent on critical elastic energy release rate of a brittle material, G c . Reader can refer to [33] for further information about calculation of s c using G c . During the PD analysis, the stretch, , between pairs of material points x i and x j are monitored in each time iteration such that if it exceeds its critical value, then the bond is failed, i.e., µ = 0, and not reevaluated in the next time steps. A demonstration of µ function is presented in Figure 2c. For a given material point's horizon, the contributions of µ lead to the φ parameter take a value between 0 and 1. If the φ = 0, then all bonds are broken, otherwise at least there is an intact bond in the horizon. A demonstration of crack integration to asymmetric distribution of particles in a PD horizon is illustrated in Figure 2. Initially, the φ is 0 since there is no broken bonds in Figure 2a. Then, as the half of the bonds are broken in Figure 2b, the φ value is obtained as 0.5. Overall, as a result of this vital feature of bonds to model damage mechanism, peridynamics can be used more effectively to handle topology optimization of cracked structures as compared to local continuum mechanics. To solve Equation (2) numerically, a mesh-less method (i.e., free from element connectivity) is used for spatial part of equation, as such the domain is divided into finite number of particles having an equivalent volume. Accordingly, Equation (2) is re-written in the following discrete form for a material point i as: where u i and u j are position vectors corresponding to i th and j th particles, N f is number of family members available in the horizon of each particle i, and V j is volume material point of j. This study will use quasi-static analysis so the density part of Equation (2) is neglected and by following [52] the linearized equation through domain β will be derived as : where K, d and x are global stiffness matrix, displacement vector, and force vector of all material points. K is symmetric matrix that includes all the stiffness contributions of the particles in their horizon. For explicit mathematical definition, reader can refer to recent study of Kefal et al. [49].

Optimization Statement, Optimality Criteria Method
TO problems in PD are focused on minimization of an objective function G with constraint of F as: where k is design variable varying between 0 (for void) and 1 (for the domain). During the topology optimization process, Young modulus of particle i can be represented via power law integration of design variable k as: where E v and E s are young modulus of void and solid, respectively. In addition, the p factor is a penalty parameter defined in the range of 1 ≤ p ≤ 5. Besides, k i value can vary from minimum design variable k min up to unity. For the current study, k min is considered as 10 −3 whereas the penalty parameter is set to as p = 3.
Compliance is an objective function that is used generally in TO with a targeted volume fraction. Here, compliance can be defined from Equation (12), thus, for the general optimization scheme, the constraints are as follows: where minimization method is subjected to the global displacements solution, updated strain field and maximum allowable volume constraintV, i.e., ratio of the target volume to the total volume of the design space.
TO method adopted in the present numerical simulations is the OC method. Essentially, this approach updates the design variable at each iteration via the following criteria: where move limit is defined as mv = 0.2 and optimality condition B i can be written as: where λ is Lagrangian multiplier chosen to satisfy the volume constraint. The sensitivity of material volume is equal to 1, and sensitivity of objective function To avoid instabilities and checkerboard patterns [53], a smooth filtering scheme is utilized during each iteration step of the OC-based TO process. Based on the family of particle i, the sensitivities are updated using the Shepard interpolation scheme defined as: where ∂Û ∂k i is updated sensitivity, r ij is distance between particle i and j, r min is filter radius, n is number of particles and ψ(r ij ) is weighting function. Here r min is set to 3.
At each iteration step, convergence of OC is monitored, and the satisfaction of the target volume constraint is insured as well. Based on error toleranceτ, the convergence criteria is defined and compared at each step as: where C m represent objective function at iteration m.

Results and Discussion
This section investigates effects of crack/damage integration for designing lightweight structures using topology optimization. Firstly, structures with no cracks are studied using both PD and FEM based topology optimization algorithms and the obtained results are compared with each other. After introducing a crack in the design domain for each case study, the TO simulations are performed again to compare the capability of PD and FEM in terms of strain energy minimization of cracked structures. For comparison of strain energy predictions, the final topology obtained from local and non-local approaches are analyzed utilizing a commercial FEM software, COMSOL, where both displacements and total strain energies of each case are calculated. During this simulation in COMSOL, the young modulus of void particles/finite elements is considered using the following penalization method: where E v and E s are young modulus of the void and solid elements. Here, elastic modulus of solid particles/elements is always considered to be 200 GPa. The penalization parameter is equal to p = 3 and Poisson's ratio is ν = 1/3. The optimization process begins with the total solid particles in the design domain and ends after satisfying the proposed convergence criterion given by Equation (20). For each model, the constraint boundary conditions, the loading condition and number of particles/elements are presented. A convergence study is performed to identify best number of particles/elements in PD/FEM. For sake of clarity in convergence results, one can refer to Appendix A where the convergence of topology optimization with respect to particles/elements sizing is provided for sample III. In addition, Section 3 includes the effect of particle sizing on the topology optimization of cracked structures. Crack/damage in PD are defined by breaking bonds according to Equations (9) and (10). Similarly for FEM, when a crack is considered, the elements corresponding to the crack region are removed from the design domain and the TO is performed accordingly.
Overall, each subsection includes investigation on structures with/without cracks. Sections 3.1 and 3.2 investigate a simply supported beam and a cantilever beam. By studying these structures, accuracy and capability of PD compared to FEM are demonstrated. In addition, capability of PD to provide better support even in a simple geometry for cracked structure is presented. Furthermore, a geometrically more complex structure such as a bracket topology is examined in Section 3.3, where divergence of FEM from PD especially in the cracked case is observed. The last benchmark problem is a bridge structure in Section 3.4. For this problem, PD-TO is introduced as a tool to provide better support in the bridge domain by utilizing cracks to reduce the strain energy. To this end, first, stress and strain energy distribution of structure without crack is obtained. Subsequently, for finding possible crack position, hot-spots with highest stress concentrations are examined. At last, by including crack/cracks in the hot-spot regions, effect of crack integration in topology optimization is studied to decrease the overall strain energy as compared to the pristine structure.

Sample I-A Simply Supported Beam
First model is a simple benchmark model as depicted in Figure 3. The structure is a simply supported beam with ratio 1 : 2 ( L = 1 m ). Namely, it is fixed from one corner and simply supported on other corner on the down edge as shown in First, in the structure without cracks, maximum vertical displacement and total strain energy obtained from PD and FEM are compared in Table 1 and Figure 4. As it can be seen from this comparison, there is a very good match between PD and FEM results with percent error of maximum displacement and total strain energy less than 0.46%. Also, the final topologies and displacement contours produced by PD are almost identıcal to those of FEM, thereby verifying the high predictive capability of PD for topology optimization of structures without crack. Later for the same case, cracked structure is investigated by keeping same volume fraction, boundary constraints and loading conditions. Since the domain itself has symmetry in geometry and constraint boundary condition, the crack is located in middle of bottom edge with symmetrical distance from corners of this edge. As can be seen from Figure 5, there is a small difference between the final topologies generated by two methods. It worth to notice that the optimum topologies are fully symmetric as the crack and geometry features possess symmetry. For a quantitative comparison of the results, the total strain energy and the maximum vertical displacement are listed in Table 2. Here, the percent difference between FEM and PD for displacement and energy is obtained as 4.5 and 5.3%, respectively. In addition, it can be observed that FEM based TO produces an optimum topology that does not support the crack tip as much as what PD-TO does ( Figure 5). Despite being a simple geometry for TO, when cracks are included, the FEM and PD optimal topology diverge from each other. The comparison of strain energies also reflects this fact such that FEM-TO results yields a larger total strain energy than PD-TO (Table 2). Hence, the better capability of PD for TO of cracked structures is revealed. Such minimization of mechanical response can also be observed for the maximum displacement of two topologies. For instance, according to the contour distribution in Figure 5, PD-TO and FEM-TO produces similar total displacement distribution, however the maximum displacement obtained by PD-TO is approximately 4.5% lower than of FEM (Table 2. Hence, this case study proves that regardless of geometrical simplicity, in case a crack is located in the structure, PD-TO produces stiffer geometry than that of FEM-TO.

Sample II-A Cantilever Beam
This section investigates another benchmark structure, which is a cantilever beam shown in Figure 6. This beam is fully clamped on left edge and a force F = 100 N is applied at the center of free edge on the right side. For PD/FEM models, the discretization size along x and y direction are 80 and 40, in the given order. The crack is considered to be in middle of top edge and with height of d = 0.2 m for cracked structure. For topology optimization of structures without crack, this value is set to 0. Here, the volume ratioV is considered to be 0.5 for both structures with/without cracks. Similar to Section 3.1, for this simple geometry, the structure without crack is used for validation, later by including a crack in design domain, capability of PD is compared to FEM L 2L F crack d L/2 Figure 6. Sample II-a cantilever beam, boundary condition and geometry. Figure 7 shows the optimum topologies obtained by FEM and PD, and associated vertical displacement contours. For a comprehensive comparison, the maximum displacement and total strain energy are listed in Table 3. As it can be seen, there is a very close agreement between both methods and the error is less than 0.66%. In the structure without cracks, the obtained geometries are symmetry since the strain energy are always distributed symmetrically during optimization steps.  Similar to previous case study, in PD, the crack is located regarding Equation (9) and in FEM, elements are excluded from the damaged/cracked region. Figure 8 shows comparison of final topologies obtained from both methods and corresponding displacement contours. As it can be seen, there is a slight difference between final topologies and the PD again provides better support for damaged area. In addition, it produces lower displacement and strain energy as well. More quantitative comparison can be observed from Table 4, where the maximum displacement and total strain energy difference between two methods is calculated as 2.9 to 3.9% respectively. On the other hand, FEM-TO only deals with constraints of the design domain without considering crack properties. Since PD include crack directly into the governing equations, it can handle cracked structures more effectively. Essentially, FEM-TO topology diverges from PD-TO one if a crack is included even in a simple geometry. It can be indicated that PD-TO produces an optimum topology having lower strain energy and displacement. In addition, PD-TO optimized topology includes better support in the cracked regions, thereby providing a less compliant structure for crack propagation as compared to FEM.

Sample III-A Bracket
In this section, geometrical complexity is enhanced to investigate the PD/FEM performance on topology optimization of cracked structures. Basically, as shown in Figure 9, an L-shaped bracket with size of (2 L = 2 m) is chosen for the topology optimization process. The bracket is fully fixed at its top edge, and a concentrated force of F = 100 N is applied at the top corner of right edge. This bracket has a notched area in its design domain which produces high concentrated hot-spot stress regions during structural analysis. The crack is included near notched area with the length of d = 0.2 m. For structure without crack, this value is set to 0. To generate PD/FEM model, the particle spacing or element size in both directions of the design domain is set to L/30.
Here, the volume ratioV is 0.5 for both structures with/without cracks. Similar to previous sections, structure without crack is used for verification purposes and the cracked case is utilized to compare the capability of PD-TO versus FEM-To for handling TO problems with cracks. In addition, the effect of discretization size on the optimum topology is discussed by including smaller mesh/element sizing i.e., L/50. For the validation, Figure 10 shows the optimum topology obtained from FEM and PD, and related displacement/strain energy. A thorough comparison for total strain energy and maximum vertical displacement between FEM-TO and PD-TO are presented in Table 5. From this comparison of figures and values in table, it can be indicated that, PD can estimate an optimum topology that is almost indistinguishable from that of FEM. In particular, the percent difference for displacement and energy between two models is less than 0.75%. Hence, the robustness of PD-TO is validated for a complex geometry without crack.  First cracked structure is discretized with 2700 particles/elements, loading and boundary condition are kept similar to those of the structure without crack. Figure 11 shows a comparison between topologies obtained from these two methods as well as their corresponding displacements. As it can be seen, when the structure become a little bit complex, the obtained topologies from PD and FEM have significant differences. In addition, displacement contours are following same trend and maximum displacement in PD is lower than the FEM one. Table 6 shows the comparison of maximum displacement and total strain energy between FEM and PD. The difference is approximately 24%, and PD has lower total strain energy and maximum displacement. Since PD-TO provides better support, it produces tougher structure. Therefore, the displacement in PD is smaller compared to FEM for the same loading condition. This case study demonstrates the advantages of PD over FEM for dealing with topology optimization of cracked structures. In fact, PD-TO provides a superior supports in order to avoid high concentrated stress regions. However, FEM-TO only satisfies the domain constraints regardless of crack location and its effective properties.  In the second cracked structure, the particles/elements number is increased to 7500 to make the model with smaller meshes to investigate effect of sizing. Figure 12 shows a comparison of topologies obtained from PD and FEM, and their corresponding displacement. Accordingly, as it can be seen, even with smaller mesh size, the FEM still diverges from PD for producing an optimum topology whereas PD provides thicker and rigid support at the vicinity of cracked regions. Table 7 provides a clear comparison of maximum displacement and total strain energy between PD and FEM models, where the percent differences are obtained approximately 18%. Since PD provides a better support in the cracked region, it has lower maximum displacement and total strain energy as compared to FEM. When the size of elements is decreased to a level that is enough to model crack as void in the FEM mesh, the FEM-TO still produces much softer topology than as what PD-TO offers. Overall, it can be concluded that PD-TO can rigorously handle topology optimization of complex cracked structures whereas FEM can not overcome this problem even using smaller meshes.

Sample IV-Bridge
Up to this section, topology optimization of simple and complex geometries are examined to reveal the advantages of peridynamics for finding an optimum topology of cracked structures. In addition, this section is devoted to establish a computational strategy for finding suitable crack locations to reduce the stress level at the hot-spots regions of the original design domain during topology optimization process. Therefore, the total strain energy of optimum topology with cracks may be smaller than the one without cracks.
A bridge like structure with length of 3L = 3 m is chosen to demonstrate the application of the computational strategy proposed herein. As depicted in Figure 13, the left and right edges of bridge are fully clamped and a uniform pressure with magnitude of F = 100 N/m is applied on the top edge of bridge. Note that these boundary conditions results in hot-spots stress regions at the constrained edges of the domain. The discretization size along x and y directions are equal to 3L/120 in both PD and FEM models. The volume fractionV is considered as 0.5 for both design domains with/without cracks. The position and size of cracks, d and w, demonstrated in Figure 13 are carefully designed through conducting an 'a priori' analysis with the known loading conditions. Here, due to the symmetry of geometry and loading/constraint conditions of bridge, the designed cracks are symmetrically integrated to the design domain. Note that if the size of crack is set to (d = 0), the structure becomes crack-free. First, structure without crack is studied to generate a reference solution for the verification purposes. Figure 14 shows the topologies and corresponding displacements contours obtained by PD-TO and FEM-TO. As seen from this figure, the PD-TO generates very similar structural geometry to the one predicted by FEM-TO with almost indistinguishable displacement contours. Moreover, the maximum displacement and total strain energy are listed and compared in Table 8. From these values, it can be indicated that the PD-TO results agree very well with those of FEM-TO with the percent error less than 1%. As validated herein, the optimum topology simply stands for a reference solution in the comparative study of structure with designed cracks.   After embedding the designed cracks in the structure, topology optimization is performed utilizing the peridynamics simulation, which generates the results shown in Figure 16a,b. When this figure is compared with Figure 14, it can be implied that the PD-TO analysis of cracked structure provides an optimum topology which supports the damaged/cracked regions better than the intact one. Also, note that the optimum topologies of both structures with/without cracks a have symmetric shape as a result of the symmetry in the strain energy distribution of bridge during TO steps. In Figure 16c,d, the stress and strain energy nearby the cracked regions are redistributed through the domain and PD-TO finally reduces total strain energy as compared to that of pristine topology. This reduction in compliance is about 3.5% as listed in Table 9. Moreover, the improve-ment in the optimum topology can be clearly recognized by comparing the maximum displacement results of the bridge with and without crack, which is about 3.3%, as listed in Table 9. Overall, it can be concluded that PD can provide an increased stiffness for topology optimization of a given structure by integrating an intended crack (i.e., having designed position and size) in the design domain. Hence, consecutive case studies presented in this study demonstrates the practical utility of PD-TO algorithms for topology optimization of cracked structures as well as designing lightweight structures through enhancing the hotspot regions via including cracks (with designed position and size) during the optimization process.

Conclusions
In this research effort, the effect of crack inclusion is studied for designing lightweight structures utilizing topology optimization method. Peridynamic theory is adopted as the main algorithm to accurately model cracks within the structure during structural analysis and optimality criteria method is employed for topology optimization process. Four different case studies with symmetrical/asymmetrical geometry and boundary conditions are solved by using PD-TO, and strain energy and displacement results are compared with those predicted by FEM-TO approach. It is demonstrated that the results of PD-TO have a very good agreement with the FEM-TO results for structures without cracks with almost negligible percent error.
To reveal the advantages of PD as compared to FEM, the same TO problems are also examined by embedding a crack with symmetric position/size in the structural domain. It is found out that the PD provides much better support for cracked structures such that total strain energy and maximum displacements estimated by PD-TO are much lower than the associated FEM results. For simple geometries with crack, the discrepancy between two methods ranges from 3 to 5%. On the other hand, this percent error increases up to 23% for a complicated structure with a crack such as a notched bracket. Moreover, this difference is not affected even if the discretization size is decreased, thereby indicating incompetence of FEM approach for handling structures with crack regardless of mesh size. Furthermore, the results show that existing initial cracks affect the final design of structures in topology optimization process, and PD-TO is much better tool for integrating crack and dealing with discontinuity like damages and internal features as compared to FEM. Overall, the main advantages of PD-TO are twofold: (i) producing optimum topologies with lower strain energy as compared to that of classical FEM-TO approach and (ii) providing better/stiffer supports for the cracked regions during the topology optimization.
Finally, a possible crack location is determined for the last case study by inspecting strain/stress distribution of the structure without crack. Then, PD-TO analysis is performed by integrating the intended crack within the structural domain. It is shown that the strain energy in the cracked structure can be reduced by proper allocation of cracks as such being even less than the original structure without cracks. Hence, it can be concluded that PD-TO can serve as a viable topology optimization tool for designing lightweight structures through wisely embedding cracks in structural domain to lower the strain energy density. Further studies can be performed to investigate the capability of PD-TO for 3D structures with surface cracks. Also, PD can be employed in structural/geometrical optimization with cracks by utilizing other optimization algorithms such as genetic algorithm in future studies.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Convergence Study for Sample III
This appendix shows the convergence study for sample III as an example for defining the best number of particles/elements for PD/FEM. By increasing number of particles from 1200 to 7500, Figure A1 shows different topology obtained from PD and FEM analysis for structures without crack. As can be seen from this figure, number of particles have a relatively small effect on the optimum topology until 2700 particles, after which shape optimization yields to the same topology. A brief comparison of maximum displacement and total strain energy convergence behaviour of both methods are presented Table A1  and Table A2. According to listed results, neither the refinement of particles nor the method being used has a very considerable influence in convergence rate of maximum displacement/total strain energy. Hence, it can be indicated that 3*30*30 discretization can be suitably used to embed cracks for both PD-TO and FEM-TO simulations, and 3*50*50 one can be utilized to check the effect discretization size on topology otpimzıation of cracked structures.