topology optimization algorithms – Global Homework Experts

Adnan Kefal1, Abdolrasoul Sohouli2, Erkan Oterkus3, Mehmet Yildiz4,5,6,*, and Afzal Suleman2
1 Faculty of Naval Architecture and Ocean Engineering, Istanbul Technical University, Maslak-Sariyer
34469, Istanbul, Turkey.
2 Department of Mechanical Engineering, University of Victoria, British Columbia, Canada
3 Integrated Manufacturing Technologies Research and Application Center, Sabanci University, Tuzla,
Istanbul 34956, Turkey
4 Department of Naval Architecture, Ocean and Marine Engineering, University of Strathclyde, Glasgow
G4 0LZ, United Kingdom
5 Composite Technologies Center of Excellence, Istanbul Technology Development Zone, Sabanci
University-Kordsa Global, Pendik, Istanbul 34906, Turkey
6 Faculty of Engineering and Natural Sciences, Sabanci University, Tuzla, Istanbul 34956, Turkey
* Corresponding author: [email protected]
Abstract: Finite element method (FEM) is commonly used with topology optimization algorithms to
determine optimum topology of load bearing structures. However, it may possess various difficulties and
limitations for handling the problems with moving boundaries, large deformations, and cracks/damages.
To remove limitations of the mesh-based topology optimization, this study presents a robust and accurate
approach based on the innovative coupling of Peridynamics (PD) (a meshless method) and topology
optimization (TO), abbreviated as PD-TO. The minimization of compliance, i.e., strain energy, is chosen
as the objective function subjected to the volume constraint. The design variable is the relative density
defined at each particle employing bi-directional evolutionary optimization approach. A filtering scheme
is also adopted to avoid the checkerboard issue and maintain the optimization stability. To present the
capability, efficiency and accuracy of the PD-TO approach, various challenging optimization problems
with and without defects (cracks) are solved under different boundary conditions. The results are
extensively compared and validated with those obtained by element free Galerkin method and FEM. The
main advantage of the PD-TO methodology is its ability to handle topology optimization problems of
cracked structures without requiring complex treatments for mesh connectivity. Hence, it can be an
alternative and powerful tool in finding optimal topologies that can circumvent crack propagation and
growth in two and three dimensional structures.
Keywords: Peridynamics, topology optimization, bi-evolutionary structural optimization, cracked

1 Introduction
The increasing need for cost-effective, light-weight, and high-performance structures have led to the
structural optimization methods classified into sizing, shape, and topology optimization [1, 2]. Topology
optimization finds an optimal material layout within a predefined design domain by minimizing or
maximizing a given objective function while satisfying design constraints. Since topology optimization
provides the most practical and potential design space among the other structural optimization methods,
it has been extensively applied in a wide range of engineering disciplines [3-7] including automotive [8-
10] and aerospace [11-14] industries. For example, various topology optimization models have been
proposed in the last few decades [15]. Remarkably, topology optimization methods are even utilized to
predict optimal distribution of shear pivot stiffness of pantographic metamaterials [16].
The homogenization method is the first approach proposed by Bendsoe [17]. Since then, several
topology optimization approaches have been introduced [18-28]. Among these, the most well-known ones
are the solid isotropic material with penalization (SIMP) method [18, 19], level set-based method [20, 21],
evolutionary structural optimization (ESO) method [22, 23], and bi-evolutionary structural optimization
(BESO) method [24-26]. Apart from these approaches, the multi-material topology optimization can also
be defined using a new definition of relative densities for multi materials as presented in [27, 28].
Additionally, the level-set method was also used for topology optimization of multi-material-based
flexoelectric composites [28]. The SIMP and BESO methods are the most commonly used approaches
among material distribution models because of their simplicity in terms of coding implementation. For
instance, the SIMP method is recently utilized to perform concurrent optimization of material density and
anisotropy in two-dimensional structures [29]. Another recent application of the SIMP approach also
involves the thermodynamic topology optimization and its comparison with heuristic optimization
schemes [30].
In the SIMP method, the design variables are defined by the relative densities, changing continuously
in the range of 1 and a nearly zero positive value, which represent the solids and voids, respectively. On
the other hand, in the ESO and BESO methods, the design variables are directly set to binary values of 1
or 0, corresponding to solid and void materials, respectively. Hence, during the optimization process,
unlike those used in the SIMP method, there is no any intermediate design variable in ESO and BESO
methods. Consequently, the ESO and BESO methods are referred to as a class of discrete topology
optimization algorithms. The ESO method is an empirical approach that gradually removes inefficient
parts from the structure during the evolution of the structure. The rejection criterion of ESO approach was

originally defined based on von Mises stress and then modified to the strain energy in [31]. The BESO
method is an extended version of the ESO approach, where the capability of the ESO method is improved
by allowing both addition of efficient parts and removal of inefficient parts simultaneously [32, 33].
Most of these aforementioned optimization techniques mainly employed a mesh-based numerical
method, i.e., finite element method (FEM) or isogeometric analysis (IGA) [34]. However, various
limitations are encountered when performing the topology optimization analysis using FEM or IGA
methods especially for problems involving moving boundaries, large deformations, and cracks/defects.
Generally, complex boundary and interface tracking algorithms and re-meshing strategies are employed
to tackle with such challenges. However, regenerating the mesh of design domain and interface tracking
in each optimization iteration can be very complex and computationally expensive [35].
Recently, meshless methods have received a great deal of attention to overcome the difficulties
mentioned in the above paragraph. In the mesh-free methods, a set of particles can arbitrarily be distributed
within the design domain without requiring any mesh connectivity, thus providing an easiness to discretize
the design domain. Nevertheless, mesh-free methods have been rarely used with the aim of the topology
optimization. To give some examples, element free Galerkin method (EFG) is one of the first meshless
methods that was employed to perform topology optimization [36-39]. Yang et al. [36] employed the EFG
method to minimize the weight of structures subjected to the displacement constraints using the SIMP
based topology optimization approach. Moreover, Shobeiri [37] also integrated the EFG method with the
BESO approach for the topology optimization of continuum structures. Shortly after, Shobeiri [38]
investigated cracked structures based on the integration of the EFG and BESO.
Radial point interpolation method [40] and local Petrov–Galerkin mixed collocation method [41] are
among other meshless approaches used in topology optimization studies. Most recently, Lin et al. [42]
proposed a structural topology optimization method using smoothed particle hydrodynamics method and
a modified SIMP based approach. Peridynamics (PD) theory is another fast-growing meshless approach
that was introduced by Silling [43] and Silling et al. [44]. The PD is considered as a nonlocal reformulation
of the classical continuum mechanics (CCM) equations [45-48]. PD has been also used for describing
growth phenomena in bones in recent years [49-51], where nonlocal biological interaction has the same
mathematical form as mechanical nonlocal interactions. The theory has been successfully used to simulate
crack growth [52], failure in multi-scale models including nanostructures [53], study damage and failure
in concrete columns [54], model fracture in functionally graded materials [55], fatigue [56], torsion [57],
buckling [58] and corrosion [59] failures, simulate hydraulic fracturing process “fracking” [60], and

predict failure in composite materials [61, 62]. PD improves the simulation of the crack growth compared
to the weakly nonlocal approaches reported in the literature [63-66]. Most recently, PD is also utilized for
toughness enhancement and material design of brittle materials by introducing micro-cracks in the vicinity
of the macro-cracks [67].
A discrete model for crack growth in continua based on the analogy with robot swarms has been
proposed recently [68-72]. In this discrete approach, the ideas of peridynamics are exploited and, when
its continuum counterpart is considered, it is proven that second and higher gradient continua are a kind
of intermediate step between nonlocal discrete interactions and PD. It has to be remarked that when
pantographic structure have infinitely stiff-in-extension constitutive bars [73-75] then their continuum
models must belong to the class which is studied in PD. Peridynamic formulation was improved for the
non-uniform discretizations and varying horizons in [76, 77]. The inherent ability of PD for dealing with
discontinuities (e.g., cracks) and their growth lends itself to topology optimization of structures in the
presence of such imperfections in a robust manner.
In this study, a novel coupling of PD and topology optimization, named concisely as PD-TO, is
introduced for the topology optimization and design of a continuum with and/or without cracks (defects).
To the best knowledge of the authors, there is no study available in the literature concerning the topology
optimization integrated with PD applied to problems with the presence of cracks. Therefore, the proposed
PD-TO approach in this study is of its first kind in the literature. The minimization of compliance, namely
strain energy, is chosen as the objective function under the volume constraint. The design variables are
the relative density defined at particles using BESO approach. The filtering scheme is adopted to avoid
the checkerboard issues and provide a high stability during optimization. The advantage of the PD-TO
methodology is its simplicity in terms of coding implementation and effectiveness, robustness, and high
accuracy to perform topology optimization of two- and three-dimensional problems without requiring any
re-meshing of the design domain. These capabilities are validated and demonstrated by solving various
challenging numerical examples.
In the following sections, the PD theory and its discretized version of equations are discussed in the
first instance (Section 2). Then, in Section 3, the BESO method is defined with an explanation of
sensitivity analysis and the filtering scheme. After that, in Section 4, the proposed PD-TO approach is
applied to different case studies including the structures with embedded cracks. Finally, the findings of
the research are summarized in the conclusions section (Section 5).

2 Peridynamics (PD)
In this study, the original Peridynamic formulation introduced by Silling [43] named as bond-based
PD is utilized due to its simplicity. However, it is straightforward to extend it to more advanced PD
formulations including state-based PD. The PD is a non-local continuum mechanics formulation.
Therefore, the continuum (analysis domain) is represented by the ensemble of infinitely small volumes
called material points (particles). As opposed to CCM, in PD, a material point can interact with other
points which are far from each other in a non-local manner and each interaction is called as Peridynamic
bond. The range of interactions is named as horizon, which represents the family of the material points,
(Figure 1). The horizon, H, can be a sphere or a circle with a radius denoted by (Figure 1) for a threedimensional (3-D) or two-dimensional (2-D) body, respectively.
Figure 1: The family
H of the material point i within the horizon radius.
The PD equation of motion of a material point located at
𝐱 can be expressed as [78]
 , , , ,         
x u x f u u x x b x t t dH t
𝜌 is the density, u is the acceleration, 𝐛 is the body force, and f is the pairwise force density
between material points
𝐱 and 𝐱′. The pairwise force density function can be defined as:
        

f u u x x f , , t f ξ η
ξ η
𝛏 represents the relative position vector between material points 𝐱 and 𝐱′ in the original
configuration, i.e.
𝛏 ൌ 𝐱െ 𝐱, and 𝛈 represents the relative displacement vector, i.e. 𝛈 ൌ 𝐮െ 𝐮. The
magnitude of the force density,
𝑓, can be expressed for an elastic isotropic material as:
f c s (3)
𝑐 is the bond constant and 𝑠 is the stretch. The bond constant for an isotropic material can be
expressed in terms of the material constants of CCM as:
 

𝐸 and representing the elastic modulus and the thickness of the 2-D domain, respectively. The
stretch of a PD bond is defined as:
 
ξ η ξ
The PD force given in Eq. (3) can be obtained from the micro-potential of the PD bond which can be
written as:
        , , 1 2
w w t cs u u x x ξ (6)
ξ is the length of the PD bond. Moreover, the strain energy density of the material point located
𝐱 can be expressed in terms of micro-potential of PD bonds associated with this material point as:
  , , ,    1 2   
W t w t dH x u u x x
Finally, integrating the strain energy density over the domain of the body,
, the total strain energy of
the body can be calculated as:
 

C W t d x, (8)
Although several analytical solutions to the PD equation given in Eq. (1) are available in the literature
[79], numerical techniques, especially meshless method, are usually utilized. Therefore, to obtain the
numerical solution based on meshless method, Eq. (1) can be written in a discrete form for the material
𝑖 as:

order now
       

    
, , , ,
i i j i j i j i
x u x f u u x x b x t t V t

N f
where 𝑁
is the number of family members of the material point 𝑖, 𝑗 represents a material point inside the
horizon (family) of the material point
𝑖, and Vj denotes the volume of material point j. Similarly, the strain
energy density of the material point
𝑖 can be written in discrete form as:
   

  
, , ,
N f
i j i j i j
W t w t V x u u x x (10)
Note that, for bond-based PD, Poisson’s ratio,
, is equal to 1 3 ⁄ for 2-D isotropic materials. Elastic
modulus value,
𝐸 is specified as the harmonic average of the elastic moduli values of the two interacting
i and j as,

i j
i j
To obtain static condition, the adaptive dynamic relaxation (ADR) technique introduced in [78, 80]
has been widely used in the literature. In ADR technique, the PD equation given in Eq. (1) yields a static
solution by imposing a fictitious damping to the system. Although ADR is a powerful technique, attaining
a static condition takes a certain number of fictitious time steps which increases the computational time
and may not be very suitable for optimization problems. Therefore, direct solution approach presented by
Bobaru et al. [81] has been utilized in this study by directly equating the inertia term in Eq. (1) to 0. This
approach requires solution of a matrix equation and can provide a quicker solution in comparison to the
ADR technique.
In order to introduce the resulting final set of equations for the direct solution [57], the PD force
density given in Eq. (2) can be linearized for material point
i by considering the work done by Silling [82]
    j       i , , , , j i ij ij ij ij ij 3 ij ij ij
ij ij
t t c
ξ ξ
f u u x x f η ξ f η T η
ξ ξ
where sign
represents the dyadic product, and ξij and ηij are position vector and relative displacement
vector defined for the interaction between the particles
i and j, respectively. The matrix Tij is a
transformation matrix for associating the relative displacement with bond direction in the reference
(undeformed) configuration. Hence this matrix is defined by the direction of cosines of the bond in the
undeformed configuration as:

 

 
 
sin cos
ij ij
 

 
cos sin
ij ij

where the direction of cosines can be calculated between the material points
i and j as:

   
  
   
ij ij
ij ij
ξ ξ
Inserting Eq. (12) into Eq. (9) and setting the acceleration term to zero leads to the PD equation of motion
for a particle
i as:

 
ij ij j i
j ij
T η V b

N f
Eq. (15) can also be presented in an alternative form as:

 
       
1  
N f
ij ij j i
j ij j
T T b
ξ u
which can be written in the form of compact matrix-vector representation as:

k d b i i i (17)
where the stiffness matrix for a particle
i can be explicitly defined as:

 
i i

 
 

i iN

 
 


1 i c k T ξ (18)
with the associated displacement vector:
 
 

 

i 1 f
u u N

         
   
a f

y a
a i N
where u x y a, denotes the displacements along x and y direction for the material point i and the
material points located within its horizon.
The PD equation of motion for the quasi-static case can be obtained as a single matrix-vector form in
the global domain
by rigorously assembling the Eq. (17) for all the material points as:
KD B (21)

K k ,

D d ,

B b (22)
K , D , and B are global stiffness matrix, displacement vector, and force vector of all material
points constructed by using local stiffness matrix
ki displacement vector di and force vector bi of each
material point
i through the assembly process, . Note that the variable N denotes the total number of
particles in the physical domain,
. Boundary conditions can be imposed by introducing additional
particles around the boundaries and creating an imaginary material domain equivalent to the horizon size
[52]. In PD, external loads can be applied as a body force onto a single layer of particles, which is also
used in this study.
3 Bi-Evolutionary Structural Optimization (BESO)
3.1 Problem Statement
The topology optimization problem for the PD material domain, , may be stated as the minimization
of an objective function,
G G , with a constraint, F F , as:
 
 
 

 
  

     
0 1,
x x
is the design variable related to the density of the material point located at the position x . In the
SIMP method, the power-law interpolation scheme is defined for the elastic modulus of each particle
based on the design variable value as well as the elastic moduli of solid and void materials as:
E E E E E i i i i       v s v  p   (24)
E v and E s are the stiffness of void and solid materials, respectively, where E E s v , and i is the
design variable of particle
i . The first E v term on the right-hand-side of Eq. (24) ensures that Ei is always
positive definite. The power,
p, of the design variable is referred to as “penalization parameter” to obtain
distinct design variables.

In the BESO method, the design variable can have two states; 1 or 0, representing solid and void
materials, respectively. For void particles, enforcing
i to acquire a very small value rather than zero
enables one to set
Ev 0 in Eq. (24). This simplifies Eq. (24) to the following form:
E E i i   p s (25)

where the design variable should be selected as i i   1 or  min for solid or void particles, respectively.


Note that the minimum value of the design variable can be defined as 0 min 10 3 and the associated

penalty parameter is chosen to have a value in the range of 1 5 p . In this study, we directly used Eq.
(25) in the BESO optimization model, where the minimum design variable and the penalty parameter are

selected as
10 and p 3, respectively. The selection of small value for design variable of void

was referred to as the soft-kill BESO approach [26]. The most beneficial aspect of using a small value
instead of 0 is to avoid from removing the particles in the design domain during the optimization process.
The stiffness is the most widely used criteria in classical engineering topology optimization with a
predefined target volume. The compliance is commonly selected as the objective function, which is the
inverse of the overall stiffness of a structure. The compliance can be defined using the total strain energy
of a structure calculated in PD given by Eq. (8). Accordingly, the optimization model can be stated with
the objective to minimize this compliance,
C , with respect to each design variable, i , as:
   

  

 
 
 
  

 
i i i
i i
i i
where the minimization problem is subjected to the solution of the global displacements to update the
strain field, and the maximum allowable (target) volume constraint denoted by the symbol
V , which is the ratio of the target volume to the total volume of the design space.
3.2 Sensitivity Numbers
In the BESO method based on FEM, the elemental sensitivity number is defined as a measure to
change the material phase of element (e.g., from solid to void or vice versa) in order to minimize the
objective function. The elements can be removed or added using the relative ranking of the sensitivity of
an individual element. Applying the same procedure to PD material points with the optimization steps,
the sensitivity number can be defined based on derivative of the strain energy or the compliance with
respect to the design variable as:

 

1 i
The sensitivity number of void elements is equal to zero when the penalization coefficient tends to go to
infinity. The derivative of the compliance can be calculated analytically as:
 
 

 

i p 1 s
i i
p C
Ci is the compliance of solid particle, which can be the basis for the compliance of a particle that
may be either solid or void. This relation can be established through the multiplication of penalized design
variable and the compliance of the solid particle as:
C C i i i   p s (29)
Substituting the compliance of solid particle defined in Eq. (29),
C C i i i s   p , into the Eq. (28), and
subsequently using Eq. (28) together with Eq. (27) lead to the final form of the sensitivity number defined
     

  
i i i i 1 min
for or
where the strain energy of each particle
Ci can be calculated at each optimization iteration step after the
solution of the displacement field as:
C W V i i i x (31)
Note that this form of the
Ci has already been incorporated into the design variable of the particle i through
the definition of the elastic modulus defined by Eq. (25).

3.3 Filter Scheme and Stabilization
A filtering or smoothing technique is applied to the sensitivity numbers to eliminate the instabilities
and checkerboard patterns as in the case of SIMP approach. For filtering purpose, a circular sub-domain
with the radius of
rmin is defined at the centroid of ith particle, which includes a set of neighbour particles,
j . Neighbour particles influence the sensitivity of the ith particle according to the distance between the
rij . Using the Shepard interpolation scheme, the sensitivity number at each particle can be
defined as:
 
 
 

 
ij j
i n

i is the smoothed sensitivity number, rij ξij is the distance between the particle i and j, n is the
total number of particles in the sub-domain of
i, and rij is a weighting function defined as:
 

 

for r r
 
r r for r r

   
min min
The optimization process can be stabilized utilizing an iterative approach. Therefore, the smoothed
sensitivity numbers,
i , are arithmetically averaged using their values from the previous iteration as:

 

 

k k
k i i

where k is the current iteration number and k
is the kth averaged sensitivity numbers.
3.4 Optimization evolution
An effective algorithm is applied to add and remove the particles to satisfy the predefined target
volume based on the ranking of the sensitivity numbers in the BESO method. The volume of the structure
is initialized with the total percentage of design domain and gradually decreased until it satisfies the target
volume. The evolutionary volume ratio introduced in [83, 84] is a controlling parameter to push the
percentage of structure’s volume to the target volume percentage. Hence, the evolutionary volume ratio,
RE , defines the percentage of the decrease or increase in the structure’s volume.
In each
kth iteration, the target volume of the next iteration Vtar k1 is defined based on the evolutionary
volume ratio as:
V V R tar tar E k k 1   1 (35)
Vtar is unity for the first iteration. In this study, the design variables are updated by means of the
optimality criteria. The pseudo code of the optimality criteria is shown in Algorithm 1. The minimum and
maximum strain energy or compliance of particles are assigned to
low and high variables, given in
Algorithm 1, at each iteration. A threshold,
th , is defined according to the ranking of the sensitivity
numbers and the target volume percentage. Accordingly, the solid particles are removed when their
sensitivity numbers are higher than the threshold. Moreover, the void particles are switched to the solid
particles if their sensitivity numbers are lower than the threshold.
𝑙𝑜𝑤 ൌ minሺ𝐶ሻ
ℎ𝑖𝑔ℎ ൌ maxሺ𝐶ሻ
𝑤ℎ𝑖𝑙𝑒 ൭൬ℎ𝑖𝑔ℎ െ 𝑙𝑜𝑤 ℎ𝑖𝑔ℎ ൰ ൐ 10
ିହ൱ 𝑡ℎ𝑒𝑛
௧௛ ൌ ൬ℎ𝑖𝑔ℎ ൅ 𝑙𝑜𝑤 2 ൰
𝑓𝑜𝑟 ሺ𝑖 ൌ 1 𝑡𝑜 𝑁ሻ 𝑡ℎ𝑒𝑛
𝑖𝑓 ሺ𝐶
൐ 𝛼௧௛ሻ 𝑡ℎ𝑒𝑛
ൌ 1
𝑒𝑙𝑠𝑒 𝑖𝑓 𝑡ℎ𝑒𝑛
ൌ 0.001 𝑜𝑟 𝜅௠௜௡
𝑒𝑛𝑑 𝑖𝑓
𝑖𝑓 ൫ሺ𝑠𝑢𝑚ሺ𝜅
ሻ െ 𝑉തሻ ൐ 0 ൯𝑡ℎ𝑒𝑛
𝑙𝑜𝑤 ൌ 𝛼
𝑒𝑙𝑠𝑒 𝑖𝑓 𝑡ℎ𝑒𝑛
ℎ𝑖𝑔ℎ ൌ 𝛼
𝑒𝑛𝑑 if
end for
𝑒𝑛𝑑 𝑤ℎ𝑖𝑙𝑒
Algorithm 1: Pseudo code of the sensitivity threshold.
Finally, a convergence criterion can be defined to stop the optimization process after satisfying the
target volume. Using a predefined error tolerance
, the convergence criterion, , can be evaluated for the
change of the objective function as:

 

 

m m
m k

 

m k

 
  
k for k
where the
C m corresponds to the objective function of mth iteration. The optimization process stops when
the convergence criterion is fully satisfied meaning that the global minimum value of the objective
function as well as the final optimized topology is achieved.
4 Numerical examples
In this section, several numerical examples are considered and solved to demonstrate the suitability
of PD-TO approach for topology optimization of structures with and without cracks. Firstly, two- and
three-dimensional topology optimization of a cantilever beam is performed for validation and parameter
selection of the present approach. Then, a cantilever beam with an edge crack is analysed to investigate
the effect of crack location on the optimal topology. Next, the effect of crack size on the topology
optimization is examined through a simply supported beam with an edge crack. The same problem is
studied to demonstrate the effect of multiple cracks at the supports as well. Finally, a simply supported
beam with an internal crack is analysed to study the effect of crack orientation on the optimal topology.
Note that, for all the cases, the objective is to design the stiffest structure based on the minimal target
volume constraint. In the following examples, unless otherwise stated, the penalization parameter is equal
𝑝 ൌ 3 and Poisson’s ratio is 1 /3. Young’s modulus of void materials is defined based on the
relation given by Eq. (25). The evolutionary ratio of the volume is set to be
RE 2% . The optimization
process begins with the total solid particles in the design domain and ends after satisfying the proposed
convergence criterion given by Eq. (36).
4.1 A cantilever beam
In this section, a cantilever beam, which is frequently reported in the literature for the validation
purposes, is analysed using the PD-TO approach. As depicted in Figure 2, design space of the beam
consists of two- and three-dimensions for
l w / 2 , where the length and width are l=100 mm and w=50
mm, respectively. The thickness is
t=7 mm in 3-D model depicted in Figure 2(b). The beam is subjected
to a concentrated force of
F=100 N at the centre of the free end. Young’s modulus of solid material is
E GPa 200 and the volume fraction constraint is set to V 0.5 . The results of the proposed approach is
compared with the optimal results of the BESO approach based on the FEM reported in [26]. The 2-D
model is discretized into
60 ൈ 30 (coarse discretization), 100 ൈ 50 (medium discretization), and 200 ൈ
(fine discretization) particles. In this comparison, the filter radius is equal to rmin 3 for all three
discretizations. The topology starts from the full solid materials in the design domain and gradually
evolves based on the BESO approach. The black and white areas denote the solid and void materials,
respectively. It is clear that the design space expands by increasing the number of particles within the
beam (Figure 3). Accordingly, the optimal topology is more detailed with the better link connections in
the structure for the fine discretization (Figure 3e-f). It can be seen from Figure 3 that the Peridynamics
based topology optimization (PD-TO) approach generates results in very close agreement with the results
obtained by FEM based topology optimization (FEM-TO).
(a) (b)
Figure 2: Design domain of a cantilever beam: (a) 2-D design domain; (b) 3-D design domain.
(a) PD 60×30 (b) FEM 60×30
(c) PD 100×50 (d) FEM 100×50
(e) PD 200×100 (f) FEM 200×100
Figure 3: Comparison of design space discretization of FEM and PD for the filter radius of 3.
The computational time for each iteration is compared in Figure 4, indicating that PD-TO approach
is computationally less expensive than the FEM based one. The total time of the optimization analysis and
the average time of each iteration in this analysis are given in Table 1 for all discretization and numerical
approaches. It can be seen that the PD-TO is approximately 25% faster than the FEM-TO for fine
discretization, which can be attributed to the form of stiffness matrix constructed and its associated
solution through matrix decomposition algorithm. It should be noted that both PD-TO and FEM-TO
solutions are obtained utilizing the same computational capabilities.

Figure 4: Time comparison of topology approaches based on FEM and PD for different
Table 1: Comparison of total time and average time for each iteration of topology approaches based on
FEM and PD for different discretization.

Coarse Medium Fine
Finite Element Method
Coarse Medium Fine
Total time [min]
Average time for each
iteration [sec]
1.45 1.33 3.12
1.18 1.40 2.84
1.01 1.71 4.10
1.21 1.47 4.04

To demonstrate the capabilities of the proposed algorithms further, we have also solved topology
optimization problem in 3-D space as shown in Figure 2b. The 3-D beam is discretized into 100×50×7
particles and the volume constraint is 0.5 with the filter radius of 3. Recall that the load is applied at the
centre of the free edge as given in Figure 2b. The design domain and the optimal topology are shown in
Figure 5(a) and (b), respectively. Expectedly, the 2-D and 3-D optimized configurations can have slightly
different topological features at the intersections of the links as can be seen from the comparison of Figure
3(c) and 5(b). This minuscule difference is related to the fact that the 3-D structure has a design space
along the thickness, which can bear the load and results in a more complex strain energy distribution
throughout the 3-D domain. Nevertheless, both optimal geometries have many common features including
void region formed in vicinity of the rear and front sections of the beam. In addition, they have a similar

stiff skeleton formed by solid particles around the border of the beam model. This case illustrates that the
proposed approach can be used for 3-D topology optimization.
Figure 5: The optimal topology obtained for 3-D analysis: (a) design domain; (b) optimized
In addition to 3-D analysis, in the following, the parameters of the topology optimization and PD are
investigated. Firstly, the effect of the filter radius is analysed on the optimal topology. The design space
is influenced by the filter radius, which plays an important role in the final topology. Figure 6 shows that
the optimal topology of the fine discretization with the filter radius of 6 is almost the same as the medium
discretization with the filter radius of 3. Enlarging the filter radius from 3 to 6 for the same discretization
decreases the number of geometrical links and increases the thickness of these links. Given that the
medium discretization with the filter radius of 3 and the fine discretization with filter radius of 6 leads to

the same topology, we have used the medium discretization for modelling of structure with cracks in the
next sections, thereby obtaining computationally less expensive solutions.
Figure 6: Analysis of the filter radius for the fine discretization (200×100) and: (a) filter radius
of 3; (b) filter radius of 6.
The second important parameter investigated herein is the horizon size. In Figure 7, the optimal
topologies are shown for the different horizon sizes
x where ranges from 2 to 5 and x is the
particle spacing for medium discretization. It is apparent from Figure 7 that as the horizon size increases,
the number of links in the topology becomes unchanged. Considering the trade-off between the
computational time and the optimum topology, we have chosen horizon size of
  3 x in the rest of this
Figure 7: Optimal topologies for 100×50 material points for different horizon sizes: (a) 2; (b) 3;
4 ; (d) 5 .
Apart from filter radius and horizon size, the effect of the filtering scheme is also investigated for the
structured and unstructured discretization. Particle based methods with uniform particle distributions can
easily discretize the complex design space commonly encountered in industrial applications. However, in
order to demonstrate that the PD-TO approach is robust and independent of uniformity of particle
distribution, we have also considered unstructured configurations as well. To make a meaningful
assessment, the design space is discretized into 5000 non-uniform and
100 50 uniform particle
distributions with the same filter radius of 3. The unstructured particle distribution is generated by using
the centroids of polygonal finite element mesh generated utilizing the open source code provided by

Talischi et al. [85]. Figure 8 shows the comparison of the optimal topologies with and without using the
sensitivity filtering for the structured and unstructured discretization. The filtering scheme can remove
significantly the checkerboard phenomena in the structured discretization of the proposed mesh-free
approach. The unstructured discretization shows less sensitivity to the filtering scheme and has less
isolated particles. The optimal results of the structured and unstructured discretization are almost
analogous to each other and the filtering technique helps to stabilize the final topology for both cases.
Hence, the filtering scheme is applied to all the test cases in the rest of this study.

Figure 8: Comparison of particle discretization and filtering effects: Structured particles (a) without
filtering, (b) with filtering; Unstructured particles (c) without filtering; (b) with filtering.
4.2 A cantilever beam with an edge crack
In this example, the effect of the crack location is investigated for a cantilever beam. The beam has a
length to width ratio of 3:1 and includes a crack with the size of
a mm 6 , which is located at a position
l measured from the clamped edge as shown in Figure 9. Two different crack locations namely,
l mm 50 and l mm 5 , are studied. The beam is subjected to a concentrated force of F=100 N at the
right bottom corner. The dimension of the beam is
150 ൈ 50 𝑚𝑚which is discretized into 90 ൈ 30
particles and Young’s modulus of solid particles is E GPa s 100 . The target volume constraint is V 0.5
for the design domain. The PD-TO simulation parameters for this test case are the same as the ones
selected in the Section 4.1.
The optimal topologies of the beam are obtained for two different positions of the crack as well as for
the beam with no crack. Figure 10 shows the histories of the compliance where the material removal
causes increase in the compliance at the initial iterations. Once the material removal is complete, the
volume fraction is fixed and the compliance starts to converge rapidly towards an optimal level for each
case. The optimal topologies of the beam are shown in Figure 11 for the beam without crack and the beam
with different crack locations. Our results are in agreement with those reported in the reference [38] where
the optimization was based on the BESO approach integrated with EFG method. In line with nature of the
optimization process, the material removal take place in low strain energy locations at each iteration.
Accordingly, it can be seen from Figure 11 that the presence of crack and/or its location alters optimized
topologies due to influence of defect (crack) on the strain energy distribution (Figure 12) in the structure.
Otherwise stated, the optimal topologies are formed such that the crack does not propagate in the beam.

Overall, such problem given in this example may not be effectively analysed utilizing optimization
approaches based on classical continuum mechanics. In fact, the classical form of equation of motion for
a material point involves the divergence of stress fields, which cannot be easily calculated at the position
of discontinuity on a design domain (e.g., stress singularity at presence of crack). In the optimization
methods based on classical mechanics, this limitation may be commonly encountered when performing
topology optimization of cracked structures. We circumvent this limitation through introducing the novel
PD-TO approach. Hence, this example problem proves that the present PD-TO approach becomes a better
candidate than those of conventional counterparts, especially for finding optimal topologies of
cracked/damaged geometries.
Figure 9: The design domain for a cantilever beam with an edge crack.
Figure 10: The iteration histories of the compliance of cantilever beam with and without crack.
Figure 11: Optimal topologies for the cantilever beam: (a) structure without crack; structure with crack:
l mm 50 and (c) l mm 5 .
Figure 12: The comparison of strain energy distributions for the cantilever beam: (a) structure
without crack; structure with crack: (b)
l mm 50 and (c) l mm 5 .
4.3 A simply supported beam with an edge crack
In this case, we have modelled a simply supported beam with length to the width ratio of is 2:1 as
shown in Figure 13. A downward concentrated force with the magnitude of F=200 N is applied onto the
middle of the top free side. The material of the solid particles has a Young’s modulus of
E GPa s 200 .
The design domain is uniformly discretized into
80 ൈ 40 particles. Each particle is assigned with a design

i and the target volume constraint is chosen to be V 0.4 . The effect of the crack length in the
plate is investigated. As illustrated in Figure 13, the crack with a length of
a is embedded at the middle
of the bottom free edge. The other problem parameters, namely,
rmin 3 and   3 x , are the same as
those chosen in Section 4.1. The topologies are investigated for three different cases with varying crack
lengths and the results are compared with optimal topology obtained for the beam without any crack.
The optimization histories of the structural compliance are depicted in Figure 14, in which the
compliance increases until satisfying the volume constraint. Then, it decreases smoothly while the volume
is almost unchanged. It is clear that the increase in the crack size leads to higher compliance values. The
optimal topologies of the beams with and without crack are shown in Figure 15. The results can be

compared with the results demonstrated in the reference [38], where the current optimal results are very
analogous to their results. It can be seen that the presence of a crack significantly influences the final
optimized topology. The strain energy distributions of the optimal topologies are compared in Figure 16.
It shows that the strain energy close to the crack tips is relatively higher than other areas and in turn, the
optimal topology evolves in such a way to prevent the crack propagation. Hence, the change of the crack
length yields different optimized topologies. The speed of the convergence of the optimization process is
sufficiently high demonstrating the numerical stability of the proposed approach.
Figure 13: Design domain of a simply supported beam with an edge crack.
Figure 14: The iteration histories of the compliance of the simply supported beam with and
without edge crack.

Figure 15: Optimal topologies for the simply supported beam: (a) structure without crack; structure with
crack: (b)
a mm 9 ; (c) a mm 12 ; (d) a mm 15 .
Figure 16: The comparison of strain energy distributions for simply supported beam: (a)
structure without crack; structure with crack: (b)
a mm 9 ; (c) a mm 12 ; (d) a mm 15 .
Additionally, we have scrutinised the effect of the multiple cracks at the constraints on the optimal
topology. The same design domain (Figure 13) is studied with two interior cracks as depicted in Figure
17. The aforementioned problem setting and parameters are kept same as well. Herein, the cracks with the
length of
a mm 18 are oriented by positive and negative 45with respect to the bottom boundary edge
(Figure 17). Figure 18 indicates that the compliance converges to a constant value demonstrating that
target volume fraction is achieved during the optimization process. In Figure 19 is presented the evolution
history of the topologies for the different iteration steps where the presence of cracks at the supports leads
to a distinctively different topology in comparison to that of the structure without crack as given in Figure
15(a). In line with the previous examples, the strain energy in the optimal topology acquires higher values
close to the crack tips (Figure 20).
Figure 17: Multiple crack formations at the constraints of the simply supported beam.
Figure 18: The iteration history of the compliance for the simply supported beam with multiple

(a) (b)
(c) (d)
(e) (f)
Figure 19: Evolution of topology for the simply supported beam with multiple cracks: (a) initial design;
(b) iteration 10; (c) iteration 20; (d) iteration 30; (e) iteration 40; (f) optimal solution.

Figure 20: Strain energy distributions for Example 5.
4.4 A simply supported beam with an interior crack
Here, we have investigated the effect of crack direction on the optimal topology of a simply supported
beam through changing the crack angle,
0 ,30 ,45 ,60 ,90      , as shown in Figure 21. The problem
setting and parameters are very similar to the previous example (Section 4.3), with the differences in the
beam dimensions (
40 20 mm2 ), crack length ( a mm 14 ) and locations (l mm 10 , l mm 10.25 ),
and target volume constraint (
V 0.5 ). Figure 22 shows that crack orientation has a notable effect on the
optimal topology. Given that high strain energy values are developed in the vicinity of the crack tips, more
solid materials points are assigned to these locations rather than the interior parts of the beam as also
observed in the previous examples. Hence, the optimal topologies evolve in the direction of preventing
the crack propagation. Figure 23 presents the iteration histories of the compliance which levels off after
satisfying the target volume constraint. In addition, Figure 24 clearly reveals that the strain energies of the
optimized topologies are maximized at the locations of the loading, constrained boundaries, crack tips,
and the junction of the three main members (links) of the geometry.
The current test problem may be analysed using the optimization approaches based on finite element
method and classical continuum mechanics; however, this process would require modelling the crack as
voids in the mesh-based discretization of the design domain. For each orientation of the crack, a
cumbersome meshing and update of the discretization would be necessity, and therefore topology
optimization process may be dependent on the size of the mesh at the crack tips. Hence, the main
disadvantage of a mesh-based discretization is that defining and tracking a crack (defect) may not be as
easy as in the case of a particle-based discretization. Since the PD naturally allows the direct solution
through particle discretization and its equations apply everywhere regardless of discontinuities, the

optimization of the current geometry would be relatively easy. Hence, the results of this complex example
reveals the superior modelling, analysis, and optimization capability of PD-TO approach in comparison
to conventional FEM-TO approaches for structures with embedded cracks.
Figure 21: The design domain of the simply supported beam having an interior crack with
different orientations

Figure 22: Optimal topologies for the simply supported beam with an interior crack: (a)
   90 , 10 l mm ; (b)    60 , 10.25 l mm ; (c) 45 , 10.25   l mm; (d)
30 , 10.25   l mm ; (e) 0 , 10   l mm .
Figure 23: The iteration histories of the compliance of the simply supported beam with an
interior crack.

Figure 24: Comparison of strain energy distributions for the simply supported beam with an
interior crack: (a)
   90 , 10 l mm ; (b)    60 , 10.25 l mm ; (c)    45 , 10.25 l mm;
   30 , 10.25 l mm ; (e)    0 , 10 l mm .
5 Concluding Remarks
This study presents an alternative topology optimization methodology to find the optimal topology of
structures with and without cracks by combining Peridynamics with the BESO approach. Unlike the meshbased methods such as FEM, topology optimization of cracked structures can be easily performed utilizing
the presented PD-TO approach. This is because crack paths can be effortlessly defined in PD without any
additional treatments to maintain mesh connectivity as in the case of mesh-based methods. Moreover, a
sensitivity filtering technique is also applied to the PD-TO methodology to avoid the checkerboard
patterns and increase the stability of optimization process. Therefore, PD-TO approach can provide
sufficient numerical accuracy and stability for the purpose of topology optimization.
To demonstrate the capability and potential applicability of the proposed PD-TO approach,
challenging 2-D and 3-D benchmark test cases are solved. The results illustrated that the crack location,
orientation and length have substantial effect on the optimized topologies. In addition, it is demonstrated
that the optimal topology can always be obtained independent of simulation parameters such as the
utilization of unstructured particle distribution, particle resolution, horizon size and filter radius, among
others. Furthermore, it is shown that PD-TO approach can be computationally faster and in turn more
efficient than the FEM-based approach. Finally, the outcomes of this study shows that the current PD-TO
approach can be an effective tool in finding optimal topologies which can prevent crack propagation and
growth in engineering structures. The future directions of the current study will focus on modelling
optimization problems including large deformations and non-linear effects.
A. Kefal greatly acknowledges the financial support provided by the Scientific and Technological
Research Council of Turkey (TUBITAK) under the grant No: 218M713. A. Sohouli and A. Suleman
acknowledge the Graduate Fellowship from the NSERC Canada Research Chair and Discovery Grant

[1] T. Yamada, K. Izui, S. Nishiwaki, and A. Takezawa, “A topology optimization method based on the level
set method incorporating a fictitious interface energy,”
Computer Methods in Applied Mechanics and
vol. 199, no. 45‐48, pp. 2876‐2891, 2010, doi: 10.1016/j.cma.2010.05.013.
K. K. Choi and N.‐H. Kim, “Structural sensitivity analysis and optimization 1: linear systems,” ed: Springer

Science & Business Media, 2006.
[3] O. Sigmund, “Design of multiphysics actuators using topology optimization ‐ Part I: One‐material
Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 49‐50, pp. 6577‐
6604, 2001, doi: 10.1016/s0045‐7825(01)00251‐1.
[4] G. H. Yoon, J. S. Jensen, and O. Sigmund, “Topology optimization of acoustic‐structure interaction
problems using a mixed finite element formulation,”
International Journal for Numerical Methods in
vol. 70, no. 9, pp. 1049‐1075, May 2007, doi: 10.1002/nme.1900.
[5] A. Gersborg‐Hansen, O. Sigmund, and R. B. Haber, “Topology optimization of channel flow problems,”
Structural and Multidisciplinary Optimization, vol. 30, no. 3, pp. 181‐192, Sep 2005, doi:
[6] E. A. Kontoleontos, E. M. Papoutsis‐Kiachagias, A. S. Zymaris, D. I. Papadimitriou, and K. C. Giannakoglou,
“Adjoint‐based constrained topology optimization for viscous flows, including heat transfer,”
Engineering Optimization, vol. 45, no. 8, pp. 941‐961, Aug 2013, doi: 10.1080/0305215x.2012.717074.
[7] S. Briot and A. Goldsztejn, “Topology optimization of industrial robots: Application to a five‐bar
Mechanism and Machine Theory, vol. 120, pp. 30‐56, Feb 2018, doi:
vol. 9, no. 3‐4, pp. 245‐249, Jul 1995, doi: 10.1007/bf01743977.
[9] J. Forsberg and L. Nilsson, “Topology optimization in crashworthiness design,”
Structural and
Multidisciplinary Optimization,
vol. 33, no. 1, pp. 1‐12, Jan 2007, doi: 10.1007/s00158‐006‐0040‐z.
[10] M. Cavazzuti, A. Baldini, E. Bertocchi, D. Costi, E. Torricelli, and P. Moruzzi, “High performance
automotive chassis design: a topology optimization based approach,”
Structural and Multidisciplinary
vol. 44, no. 1, pp. 45‐56, Jul 2011, doi: 10.1007/s00158‐010‐0578‐7.
[11] D. Inoyama, B. P. Sanders, and J. J. Joo, “Topology Optimization Approach for the Determination of the
Multiple‐Configuration Morphing Wing Structure,”
Journal of Aircraft, vol. 45, no. 6, pp. 1853‐1862, Nov‐
Dec 2008, doi: 10.2514/1.29988.
[12] J. H. Zhu, W. H. Zhang, and L. Xia, “Topology Optimization in Aircraft and Aerospace Structures Design,”
Archives of Computational Methods in Engineering, vol. 23, no. 4, pp. 595‐622, Dec 2016, doi:
[13] K. Maute and M. Allen, “Conceptual design of aeroelastic structures by topology optimization,”
Structural and Multidisciplinary Optimization, vol. 27, no. 1‐2, pp. 27‐42, May 2004, doi:
[14] R. Zhang, X. Zhang, G. Lorenzini, and G. Xie, “Material combinations and parametric study of thermal and
mechanical performance of pyramidal core sandwich panels used for hypersonic aircrafts,”
Mechanics and Thermodynamics,
vol. 28, no. 6, pp. 1905‐1924, Nov 2016, doi: 10.1007/s00161‐016‐
[15] O. Sigmund and K. Maute, “Topology optimization approaches A comparative review,”
Structural and
Multidisciplinary Optimization,
vol. 48, no. 6, pp. 1031‐1055, Dec 2013, doi: 10.1007/s00158‐013‐0978‐
[16] F. dell’Isola
et al., “Advances in pantographic structures: design, manufacturing, models, experiments
and image analyses,”
Continuum Mechanics and Thermodynamics, vol. 31, no. 4, pp. 1231‐1282, Jul
2019, doi: 10.1007/s00161‐019‐00806‐x.
Computer Methods in Applied Mechanics and Engineering, vol. 71, no. 2,
pp. 197‐224, Nov 1988, doi: 10.1016/0045‐7825(88)90086‐2.
[18] M. P. Bendsoe and O. Sigmund, “Material interpolation schemes in topology optimization,”
Archive of
Applied Mechanics,
vol. 69, no. 9‐10, pp. 635‐654, Nov 1999, doi: 10.1007/s004190050248.
[19] B. M. Philip and S. Ole,
Topology optimization: Theory, methods and applications. Berlin: Springer, 2003.
[20] M. Y. Wang, X. Wang, and D. Guo, “A level set method for structural topology optimization,” Computer
methods in applied mechanics and engineering,
vol. 192, no. 1, pp. 227‐246, 2003.
[21] N. P. van Dijk, K. Maute, M. Langelaar, and F. van Keulen, “Level‐set methods for structural topology
optimization: a review,”
Structural and Multidisciplinary Optimization, vol. 48, no. 3, pp. 437‐472, Sep
2013, doi: 10.1007/s00158‐013‐0912‐y.
Computers & Structures, vol. 49, no. 5, pp. 885‐896, Dec 1993, doi: 10.1016/0045‐7949(93)90035‐c.
[23] X. D. Huang and Y. M. Xie, “A further review of ESO type methods for topology optimization,”
and Multidisciplinary Optimization,
vol. 41, no. 5, pp. 671‐683, May 2010, doi: 10.1007/s00158‐010‐
[24] O. M. Querin, V. Young, G. P. Steven, and Y. M. Xie, “Computational efficiency and validation of bi‐
directional evolutionary structural optimisation,”
Computer Methods in Applied Mechanics and
vol. 189, no. 2, pp. 559‐573, 2000, doi: 10.1016/s0045‐7825(99)00309‐6.
[25] X. Huang, Y. M. Xie, and M. C. Burry, “Advantages of bi‐directional evolutionary structural optimization
(BESO) over evolutionary structural optimization (ESO),”
Advances in Structural Engineering, vol. 10, no.
6, pp. 727‐737, Dec 2007, doi: 10.1260/136943307783571436.
[26] X. Huang and M. Xie,
Evolutionary topology optimization of continuum structures: methods and
. John Wiley & Sons, 2010.
[27] J. H. Sjolund, D. Peeters, and E. Lund, “A new thickness parameterization for Discrete Material and
Thickness Optimization,”
Structural and Multidisciplinary Optimization, vol. 58, no. 5, pp. 1885‐1897,
Nov 2018, doi: 10.1007/s00158‐018‐2093‐1.
[28] H. Ghasemi, H. S. Park, and T. Rabczuk, “A multi‐material level set‐based topology optimization of
flexoelectric composites,”
Computer Methods in Applied Mechanics and Engineering, vol. 332, pp. 47‐62,
Apr 2018, doi: 10.1016/j.cma.2017.12.005.
[29] N. Ranaivomiarana, F. X. Irisarri, D. Bettebghor, and B. Desmorat, “Concurrent optimization of material
spatial distribution and material anisotropy repartition for two‐dimensional structures,”
Mechanics and Thermodynamics,
vol. 31, no. 1, pp. 133‐146, Jan 2019, doi: 10.1007/s00161‐018‐0661‐7.
[30] D. R. Jantos, C. Riedel, K. Hackl, and P. Junker, “Comparison of thermodynamic topology optimization
with SIMP,”
Continuum Mechanics and Thermodynamics, vol. 31, no. 2, pp. 521‐548, Mar 2019, doi:
[31] P. Tanskanen, “The evolutionary structural optimization method: theoretical aspects,”
Methods in Applied Mechanics and Engineering,
vol. 191, no. 47‐48, pp. 5485‐5498, 2002, Art no. Pii
s0045‐7825(02)00464‐4, doi: 10.1016/s0045‐7825(02)00464‐4.
[32] O. M. Querin, G. P. Steven, and Y. M. Xie, “Evolutionary structural optimisation (ESO) using a
bidirectional algorithm,”
Engineering Computations, vol. 15, no. 8, pp. 1031‐+, 1998, doi:
[33] X. Y. Yang, Y. M. Xie, G. P. Steven, and O. M. Querin, “Bidirectional evolutionary method for stiffness
Aiaa Journal, vol. 37, no. 11, pp. 1483‐1488, Nov 1999, doi: 10.2514/2.626.
[34] H. Ghasemi, H. S. Park, and T. Rabczuk, “A level‐set based IGA formulation for topology optimization of
flexoelectric materials,”
Computer Methods in Applied Mechanics and Engineering, vol. 313, pp. 239‐
258, 2017.
[35] O. Sigmund and J. Petersson, “Numerical instabilities in topology optimization: A survey on procedures
dealing with checkerboards, mesh‐dependencies and local minima,”
Structural Optimization, vol. 16, no.
1, pp. 68‐75, Aug 1998, doi: 10.1007/bf01214002.
[36] X. J. Yang, J. Zheng, and S. Y. Long, “Topology optimization of continuum structures with displacement
constraints based on meshless method,”
International Journal of Mechanics and Materials in Design, vol.
13, no. 2, pp. 311‐320, Jun 2017, doi: 10.1007/s10999‐016‐9337‐2.

[37] V. Shobeiri, “Topology optimization using bi‐directional evolutionary structural optimization based on
the element‐free Galerkin method,”
Engineering Optimization, vol. 48, no. 3, pp. 380‐396, Mar 2016,
doi: 10.1080/0305215x.2015.1012076.
[38] V. Shobeiri, “The topology optimization design for cracked structures,”
Engineering Analysis with
Boundary Elements,
vol. 58, pp. 26‐38, Sep 2015, doi: 10.1016/j.enganabound.2015.03.002.
[39] Q. Z. He, Z. Kang, and Y. Q. Wang, “A topology optimization method for geometrically nonlinear
structures with meshless analysis and independent density field interpolation,”
vol. 54, no. 3, pp. 629‐644, Sep 2014, doi: 10.1007/s00466‐014‐1011‐7.
[40] J. Zheng, S. Y. Long, Y. B. Xiong, and G. Y. Li, “A Topology Optimization Design for the Continuum
Structure Based on the Meshless Numerical Technique,”
Cmes‐Computer Modeling in Engineering &
vol. 34, no. 2, pp. 137‐154, Sep 2008.
[41] S. Li and S. N. Aturi, “Topology‐optimization of structures based on the MLPG mixed collocation
Cmes‐Computer Modeling in Engineering & Sciences, vol. 26, no. 1, pp. 61‐74, Mar 2008.
[42] J. Lin, Y. Guan, G. Zhao, H. Naceur, and P. Lu, “Topology optimization of plane structures using smoothed
particle hydrodynamics method,”
International Journal for Numerical Methods in Engineering, 2016.
[43] S. A. Silling, “Reformulation of elasticity theory for discontinuities and long‐range forces,”
Journal of the
Mechanics and Physics of Solids,
vol. 48, no. 1, pp. 175‐209, Jan 2000, doi: 10.1016/s0022‐
[44] S. A. Silling, M. Epton, O. Weckner, J. Xu, and E. Askari, “Peridynamic states and constitutive modeling,”
Journal of Elasticity, vol. 88, no. 2, pp. 151‐184, Aug 2007, doi: 10.1007/s10659‐007‐9125‐1.
[45] F. Dell’Isola, U. Andreaus, and L. Placidi, “At the origins and in the vanguard of peridynamics, non‐local
and higher‐gradient continuum mechanics: An underestimated and still topical contribution of Gabrio
Mathematics and Mechanics of Solids, vol. 20, no. 8, pp. 887‐928, Sep 2015, doi:
[46] F. Dell’Isola, A. Della Corte, R. Esposito, and L. Russo, “Some Cases of Unrecognized Transmission of
Scientific Knowledge: From Antiquity to Gabrio Piola’s Peridynamics and Generalized Continuum
Generalized Continua as Models for Classical and Advanced Materials, vol. 42, pp. 77‐128,
2016, doi: 10.1007/978‐3‐319‐31721‐2_5.
[47] S. Ebrahimi, D. J. Steigmann, and K. Komvopoulos, “PERIDYNAMICS ANALYSIS OF THE NANOSCALE
Journal of Mechanics of
Materials and Structures,
vol. 10, no. 5, pp. 559‐572, Dec 2015, doi: 10.2140/jomms.2015.10.559.
[48] M. Taylor and D. J. Steigmann, “A two‐dimensional peridynamic model for thin plates,”
and Mechanics of Solids,
vol. 20, no. 8, pp. 998‐1010, Sep 2015, doi: 10.1177/1081286513512925.
[49] T. Lekszycki and F. dell’Isola, “A mixture model with evolving mass densities for describing synthesis and
resorption phenomena in bones reconstructed with bio‐resorbable materials,”
Zamm‐Zeitschrift Fur
Angewandte Mathematik Und Mechanik,
vol. 92, no. 6, pp. 426‐444, Jun 2012, doi:
[50] I. Giorgio, U. Andreaus, D. Scerrato, and F. dell’Isola, “A visco‐poroelastic model of functional adaptation
in bones reconstructed with bio‐resorbable materials,”
Biomechanics and Modeling in Mechanobiology,
vol. 15, no. 5, pp. 1325‐1343, Oct 2016, doi: 10.1007/s10237‐016‐0765‐6.
[51] I. Giorgio, U. Andreaus, F. Dell’Isola, and T. Lekszycki, “Viscous second gradient porous materials for
bones reconstructed with bio‐resorbable grafts,”
Extreme Mechanics Letters, vol. 13, pp. 141‐147, May
2017, doi: 10.1016/j.eml.2017.02.008.
[52] S. A. Silling and E. Askari, “A meshfree method based on the peridynamic model of solid mechanics,”
Computers & Structures, vol. 83, no. 17‐18, pp. 1526‐1535, Jun 2005, doi:
[53] F. Bobaru, S. A. Silling, and H. Jiang, “Peridynamic fracture and damage modeling of membranes and
nanofiber networks,” in
XI Int. Conf. Fract., Turin, Italy, 2005.
[54] F. Shen, Q. Zhang, and D. Huang, “Damage and Failure Process of Concrete Structure under Uniaxial
Compression Based on Peridynamics Modeling,”
Mathematical Problems in Engineering, 2013, Art no.
631074, doi: 10.1155/2013/631074.
[55] Z. Q. Cheng, G. F. Zhang, Y. N. Wang, and F. Bobaru, “A peridynamic model for dynamic fracture in
functionally graded materials,”
Composite Structures, vol. 133, pp. 529‐546, Dec 2015, doi:
[56] G. F. Zhang, Q. Le, A. Loghin, A. Subramaniyan, and F. Bobaru, “Validation of a peridynamic model for
fatigue cracking,”
Engineering Fracture Mechanics, vol. 162, pp. 76‐94, Aug 2016, doi:
Journal of Mechanics of Materials and Structures, vol. 10, no. 2, pp. 167‐193, Mar 2015, doi:
[58] B. Kilic and E. Madenci, “Structural stability and failure analysis using peridynamic theory,” (in English),
International Journal of Non‐Linear Mechanics, Article vol. 44, no. 8, pp. 845‐854, Oct 2009, doi:
[59] D. De Meo, C. Diyaroglu, N. Zhu, E. Oterkus, and M. A. Siddiq, “Modelling of stress‐corrosion cracking by
using peridynamics,” (in English),
International Journal of Hydrogen Energy, Article vol. 41, no. 15, pp.
6593‐6609, Apr 2016, doi: 10.1016/j.ijhydene.2016.02.154.
[60] S. Nadimi, I. Miscovic, and J. McLennan, “A 3D peridynamic simulation of hydraulic fracture process in a
heterogeneous medium,” (in English),
Journal of Petroleum Science and Engineering, Article vol. 145, pp.
444‐452, Sep 2016, doi: 10.1016/j.petrol.2016.05.032.
[61] B. Kilic, A. Agwai, and E. Madenci, “Peridynamic theory for progressive damage prediction in center‐
cracked composite laminates,”
Composite Structures, vol. 90, no. 2, pp. 141‐151, Sep 2009, doi:
[62] Y. L. Hu, N. V. De Carvalho, and E. Madenci, “Peridynamic modeling of delamination growth in composite
Composite Structures, vol. 132, pp. 610‐620, Nov 2015, doi:
[63] L. Placidi and E. Barchiesi, “Energy approach to brittle fracture in strain‐gradient modelling,”
of the Royal Society a‐Mathematical Physical and Engineering Sciences,
vol. 474, no. 2210, Feb 2018, Art
no. 20170878, doi: 10.1098/rspa.2017.0878.
Mathematics and
Mechanics of Complex Systems,
vol. 6, no. 2, pp. 77‐100, 2018, doi: 10.2140/memocs.2018.6.77.
[65] L. Placidi, A. Misra, and E. Barchiesi, “Two‐dimensional strain gradient damage modeling: a variational
Zeitschrift Fur Angewandte Mathematik Und Physik, vol. 69, no. 3, Jun 2018, Art no. 56, doi:
[66] L. Placidi, A. Misra, and E. Barchiesi, “Simulation results for damage with evolving microstructure and
growing strain gradient moduli.,” ed: Continuum Mechanics and Thermodynamics, 2018, pp. 1‐21.
[67] M. F. Basoglu, Z. Zerin, A. Kefal, and E. Oterkus, “A computational model of peridynamic theory for
deflecting behavior of crack propagation with micro‐cracks,”
Computational Materials Science, vol. 162,
pp. 33‐46, May 2019, doi: 10.1016/j.commatsci.2019.02.032.
[68] A. Della Corte, A. Battista, F. Dell’Isola, and I. Giorgio, “Modeling deformable bodies using discrete
systems with centroid‐based propagating interaction: fracture and crack evolution.,” ed. Springer,
Singapore: Mathematical Modelling in Solid Mechanics, 2017, pp. 59‐88.
[69] A. Della Corte, A. Battista, and F. dell’Isola, “Referential description of the evolution of a 2D swarm of
robots interacting with the closer neighbors: Perspectives of continuum modeling via higher gradient
International Journal of Non‐Linear Mechanics, vol. 80, pp. 209‐220, Apr 2016, doi:

[70] J. Wiech, V. A. Eremeyev, and I. Giorgio, “Virtual spring damper method for nonholonomic robotic
swarm self‐organization and leader following.,” vol. 30, ed: Continuum Mechanics and
Thermodynamics, pp. 1091‐1102.
[71] A. Battista, L. Rosa, R. dell’Erba, and L. Greco, “Numerical investigation of a particle system compared
with first and second gradient continua: Deformation and fracture phenomena*,”
Mathematics and
Mechanics of Solids,
vol. 22, no. 11, pp. 2120‐2134, Nov 2017, doi: 10.1177/1081286516657889.
[72] R. dell’Erba, “Swarm robotics and complex behaviour of continuum material.,” ed: Continuum
Mechanics and Thermodynamics, 2018, pp. 1‐26.
[73] E. Turco, F. dell’Isola, A. Cazzani, and N. L. Rizzi, “Hencky‐type discrete model for pantographic
structures: numerical comparison with second gradient continuum models,”
Zeitschrift Fur Angewandte
Mathematik Und Physik,
vol. 67, no. 4, Aug 2016, Art no. 85, doi: 10.1007/s00033‐016‐0681‐8.
[74] E. Turco, F. dell’Isola, N. L. Rizzi, R. Grygoruk, W. H. Muller, and C. Liebold, “Fiber rupture in sheared
planar pantographic sheets: Numerical and experimental evidence,”
Mechanics Research
vol. 76, pp. 86‐90, Sep 2016, doi: 10.1016/j.mechrescom.2016.07.007.
[75] J. J. Alibert, P. Seppecher, and F. Dell’Isola, “Truss modular beams with deformation energy depending
on higher displacement gradients,”
Mathematics and Mechanics of Solids, vol. 8, no. 1, pp. 51‐73, Feb
2003, doi: 10.1177/108128603029658.
[76] H. L. Ren, X. Y. Zhuang, Y. C. Cai, and T. Rabczuk, “Dual‐horizon peridynamics,”
International Journal for
Numerical Methods in Engineering,
vol. 108, no. 12, pp. 1451‐1476, Dec 2016, doi: 10.1002/nme.5257.
[77] H. L. Ren, X. Y. Zhuang, and T. Rabczuk, “Dual‐horizon peridynamics: A stable solution to varying
Computer Methods in Applied Mechanics and Engineering, vol. 318, pp. 762‐782, May 2017,
doi: 10.1016/j.cma.2016.12.031.
[78] E. Madenci and E. Oterkus,
Peridynamic theory and its applications. Springer, 2014.
[79] Y. Mikata, “Analytical solutions of peristatic and peridynamic problems for a 1D infinite rod,”
International Journal of Solids and Structures, vol. 49, no. 21, pp. 2887‐2897, Oct 2012, doi:
[80] P. Underwood, “Dynamic relaxation(in structural transient analysis),” ed: Computational methods for
transient analysis(A 84‐29160 12‐64), 1983, pp. 245‐265.
[81] F. Bobaru, M. J. Yang, L. F. Alves, S. A. Silling, E. Askari, and J. F. Xu, “Convergence, adaptive refinement,
and scaling in 1D peridynamics,”
International Journal for Numerical Methods in Engineering, vol. 77, no.
6, pp. 852‐877, Feb 2009, doi: 10.1002/nme.2439.
[82] S. A. Silling, “Linearized Theory of Peridynamic States,”
Journal of Elasticity, vol. 99, no. 1, pp. 85‐111,
Mar 2010, doi: 10.1007/s10659‐009‐9234‐0.
[83] X. Huang, Y. M. Xie, and M. C. Burry, “A new algorithm for bi‐directional evolutionary structural
optimization,” vol. 49, ed: SME International Journal Series C Mechanical Systems, Machine Elements
and Manufacturing, 2006, pp. 1091‐1099.
[84] X. Huang and Y. M. Xie, “Convergent and mesh‐independent solutions for the bi‐directional evolutionary
structural optimization method,”
Finite Elements in Analysis and Design, vol. 43, no. 14, pp. 1039‐1049,
Oct 2007, doi: 10.1016/j.finel.2007.06.006.
[85] C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes, “PolyMesher: a general‐purpose mesh
generator for polygonal elements written in Matlab,”
Structural and Multidisciplinary Optimization, vol.
45, no. 3, pp. 309‐328, Mar 2012, doi: 10.1007/s00158‐011‐0706‐z.