Structural Shape Optimization Based on Multi-Patch Weakly Singular IGABEM and Particle Swarm Optimization Algorithm in Two-Dimensional Elastostatics

Abstract

In this paper, a multi-patch weakly singular isogeometric boundary element method (WSIGABEM) for two-dimensional elastostatics is proposed. Since the method is based on the weakly singular boundary integral equation, quadrature techniques, dedicated to the weakly singular and regular integrals, are applied in the method. A new formula for the generation of collocation points is suggested to take full advantage of the multi-patch technique. The generated collocation points are essentially inside the patches without any correction. If the boundary conditions are assumed to be continuous in every patch, no collocation point lies on the discontinuous boundaries, thus simplifying the implementation. The multi-patch WSIGABEM is verified by simple examples with analytical solutions. The features of the present multi-patch WSIGABEM are investigated by comparison with the traditional IGABEM. Furthermore, the combination of the present multi-patch WSIGABEM and the particle swarm optimization algorithm results in a shape optimization method in two-dimensional elastostatics. By changing some specific control points and their weights, the shape optimizations of the fillet corner, the spanner, and the arch bridge are verified to be effective.

1. Introduction

Isogeometric Analysis (IGA) was first proposed by Hughes et al. in 2005 [1]. This method aims to bridge the gap between Computer-Aided Design (CAD) and Computer-Aided Engineering (CAE) through the use of the same geometric description. In many engineering fields, such as automotive, shipbuilding, and aerospace, mesh generation occupies a significant portion of the entire analysis process, accounting for about 80% of the total time. The conversion of CAD geometric data to finite element meshes requires extensive manual operations, and the adaptive mesh refinement needs to interact with CAD systems. Furthermore, frequent adjustments to the product’s shape during the design process necessitate the continual regeneration of meshes. Therefore, the simplification, automation, and intelligence of mesh generation are important directions in the development of computational mechanics and CAE software (https://www.sw.siemens.com/en-US/technology/computer-aided-engineering-cae/, accessed on 22 April 2024). In isogeometric analysis, CAE software directly uses the same geometrical model built by CAD for computation without the need for mesh conversion. Mesh refinement can be achieved by adding knots or increasing the order of the spline function [2,3]. Isogeometric analysis technology significantly improves the accuracy of the engineering analysis and simplifies workflows by making information exchange during geometric refinement more direct and efficient. Isogeometric analysis has been explored mainly on the single-patch technique [4,5,6,7,8,9,10,11,12,13,14], while relatively few studies have been conducted on the multi-patch technique. The advantage of the multi-patch method lies in its ability to be applied to irregularly trimmed NURBSs [15,16].

The Boundary Element Method (BEM) is a numerical method that converts the boundary value problem of domain differential equations into integral equations on the boundary, thereby performing numerical solutions only on the boundary. This method significantly reduces the dimensionality of the problem, potentially improving computational efficiency in certain cases [17,18]. Over the past few decades, the boundary element method has developed into a widely used and powerful computational tool in engineering fields and is often an alternative to the more commonly used finite element method. The BEM utilizes only geometric information from the boundary. Since the boundary contains all the geometric information of the entire entity, the boundary element method is more efficient in CAD/CAE communication than the finite element method. In the BEM, the numerical techniques for the singular integrals are key research focuses and challenges. In recent years, many scholars have investigated singular integral problems in various applied fields, as referenced in [19,20,21,22,23,24,25,26].

The application of isogeometric analysis in the boundary element method began with research by Simpson and others in the framework of two-dimensional static structural analysis of elastic mechanics [27,28]. At the same time, Gu and others applied it to solving three-dimensional potential problems [29]. The Isogeometric Boundary Element Method (IGABEM) has been widely used to solve various problems, such as linear elasticity [30,31,32,33], acoustics [34,35,36], fluid mechanics [37], fracture mechanics [38,39], and inclusion problems [40,41]. In recent years, the application of IGABEM in shape and topology optimization has also received attention [42,43,44]. For a comprehensive introduction and detailed techniques regarding IGABEM, readers are referred to the monograph by Bear [45]. Collocation points are essential for the completeness of IGABEM. Different scholars hold different opinions on the optimal collocation method, and the common approaches are the Gaussian collocation method [46], orthogonal function collocation method [47], and uniform collocation method [48], ref. [49] evaluates different collocation methods and concludes that the Greville abscissae collocation method is a simple and efficient method that produces quite accurate and stable results. Similarly, refs. [50,51] extend the development of collocation methods within the framework of IGA to multi-patch NURBS configurations.

In the community of researchers in the boundary element method, the boundary integral equations of weakly singular form or nonsingular form are of special interest. Gu et al. [29] developed a weakly singular isogeometric boundary element method for potential problems. This approach effectively improved the computational efficiency of the boundary element method for 3D potential problems since the computation of a large number of nearly singular and strongly singular integrals was not necessary. Heltai et al. [37] proposed a non-singular isogeometric boundary element method for solving the 3D Stokes flow problem, making it possible to use the standard numerical integration formulation to solve the equations. Wang et al. [32] proposed an efficient multi-patch nonsingular isogeometric boundary element method for 3D elastostatic problems. However, the strongly and weakly singular integrals were still involved in the implementation, i.e., the integrals in Equation (46) of reference [28]. A large number of Gauss points were still required in the numerical integration.

Due to IGA closely linking the design model to the analysis model, it is more suitable for shape optimization design [52]. NURBS splines have efficient shape parameters and can represent complex free-form shapes with fewer control points [53]. Since the BEM only requires discretization on the boundary, and IGA is efficient on the geometric description, IGABEM can significantly reduce both computational and storage requirements in shape optimization. IGABEM can easily achieve adaptive discretization by adjusting the knots and order of the spline function, which means that the number of refinements can be increased in regions that require higher precision, while the number of refinements can be kept low in less important regions, thus optimizing the use of computational resources. In general, IGABEM provides higher accuracy and efficiency in shape optimization, especially when dealing with complex geometries and high-dimensional problems. Currently, isogeometric analysis methods have been widely applied to a variety of engineering optimization problems, such as optimization of beam structures [54], shape optimization of vibrating membranes [55], and shape design optimization of SPH fluid–structure interactions [56]. Lian H J et al. [57,58] developed a gradient-based optimization algorithm that provides an efficient solution for shape optimization of 2D and 3D linear elastic structures; however, the process requires complex sensitivity analysis. Sun et al. [59] proposed a structural shape optimization method combining IGABEM and the Particle Swarm Optimization (PSO) algorithm. By using NURBS control points as design parameters, the unified description of the design model, analysis model, and optimization model was realized and avoids complicated sensitivity analysis. However, the paper does not make full use of the characteristics of NURBS splines, such as weights that affect the shape of curves are not included in the design parameters.

In this paper, a multi-patch weakly singular isogeometric boundary element method (WSIGABEM) for two-dimensional elasticity problems is proposed, aiming to provide a numerical method with higher accuracy and lower computational cost. The multi-patch WSIGABEM simplifies the numerical integration due to the absence of the strong singular integral, with only weakly singular integrals remaining. The numerical integration is evaluated by Telles [60]. A simple formula is suggested for the generation of the collocation points, which automatically locate inside the patch without any correction. By using the multi-patch technique so that the displacements and tractions in each patch are continuous, all the collocation points are continuous, i.e., there are no collocation points at any edge or corner. The numerical computation results of the multi-patch WSIGABEM are verified for accuracy by comparing with the traditional IGABEM via two simple problems with analytical solutions. Furthermore, based on the multi-patch WSIGABEM and the PSO algorithm, a shape optimization method suitable for complicated models is proposed. Both the coordinates and weights of the control points serve as the design parameters in the three shape optimization examples, which show the elegance of the present method.

2. B-Splines and NURBS

B-spline curves are generated from the control points and B-spline basis functions. Among the many different methods for the evaluation of B-spline functions, the recursive method is the most efficient in computer implementations. One can define a monotonically non-decreasing sequence of real numbers as the knot vector

$$\Xi =\left\{{\xi }_1,{\xi }_2,{\xi }_3,\dots ,{\xi }_{n+p+1}\right\}, {\xi }_i\le {\xi }_{i+1}, i=1,2,\dots n+p$$

(1)

where ${\xi }_i$ is a knot, p is the order of the basis function, and n is the number of B-spline basis functions generated by the knot vector. A knot vector is uniform if all internal knots are equidistantly distributed, otherwise $\Xi$ is non-uniform. Knot vectors with knots repeating p + 1 times at the beginning and the end are open knot vectors, which are commonly used in industrial design software. The B-spline basis function with order p > 0 is defined recursively by

$$N_{i,p}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}(\xi )+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}(\xi )$$

(2)

with

$$N_{i,0}(\xi )=\left\{\begin{array}{l}1, {\xi }_i<\xi <{\xi }_{i+1}\\ 0, others\end{array}\right.$$

(3)

A curve in geometric space can be expressed as a linear combination of the B-spline basis function and the control points, which are the selected points in the space. A B-spline curve of order p is defined as follows

$$\mathit{C}(\xi )={\displaystyle \sum _{i=1}^n{N}_{i,p}(\xi ){\mathit{P}}_i ,{\xi }_1\le \xi \le {\xi }_{n+p+1}}$$

(4)

where P_i represents the ith control point. The order represents the complexity of the curve. From Equation (2), it can be seen that the high-order basis functions are linear combinations of the basis functions of one order lower.

In the subsequent implementation of the boundary element method, the calculation of the derivatives of the basis functions, used for the approximation of the unknown field, is indispensable. The first derivative of the ith B-spline basis function of order p is defined as

$$\frac{d{N}_{i,p}(\xi )}{d\xi }=\frac{p}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}(\xi )+\frac{p}{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}(\xi )$$

(5)

High-order derivatives are obtained by the recurrence relation:

$$N_{i,p}^{(k)}=p\left(\frac{N_{i,p-1}^{(k-1)}}{{\xi }_{i+p}-{\xi }_i}-\frac{N_{i+1,p-1}^{(k-1)}}{{\xi }_{i+p+1}-{\xi }_{i+1}}\right)$$

(6)

Based on B-spline curves, NURBS curves introduce a weight w, which reflects the influence of control points on the NURBS curve. Generally, the larger the weight, the closer the control point to the curve. The NURBS curve is defined as follows

$$\mathit{C}(\xi )={\displaystyle \sum _{i=1}^n{R}_{i,p}(\xi )}{\mathit{P}}_i$$

(7)

where the NURBS basis function is defined by

$$R_{i,p}(\xi )=\frac{N_{i,p}(\xi )w_i}{{\displaystyle \sum _{j=1}^n{N}_{j,p}(\xi )w_j}}$$

(8)

The introduction of weights makes the expression of curve shapes more flexible, allowing the modification of curve shapes not only by changing the positions of control points but also by adjusting the weight.

3. Multi-Patch WSIGABEM

3.1. Weakly Singular Boundary Integral and Discretization

In a physical domain Ω with boundary Γ, the boundary integral equation of linear elastostatics in the absence of body forces is expressed as

$$C_{ij}(\mathit{\xi })u_j(\mathit{\xi })+{\displaystyle {\int }_{\Gamma }{t}_{ij}^{\ast }(\mathit{\xi },\mathit{x})u_j(\mathit{x})d\Gamma (\mathit{x})}={\displaystyle {\int }_{\Gamma }{u}_{ij}^{\ast }(\mathit{\xi },\mathit{x})t_j(\mathit{x})d\Gamma (\mathit{x})}$$

(9)

where the integral on the left side represents a Cauchy principal value integral, C_ij represents the jump term, and $t_{ij}^{\ast }(\mathit{\xi },\mathit{x})$ and $u_{ij}^{\ast }(\mathit{\xi },\mathit{x})$ are, respectively, the plane-strain fundamental solution kernels for 2D problems, i.e.,

$$t_{km}^{\ast }(\mathit{x},\mathit{\xi })=-\frac{1}{4\pi G(1-\nu )r}\left\{\frac{\partial r}{\partial n}[(1-2\upsilon ){\delta }_{km}\right.\left.+2r_{,k}{r}_{,m}]-(1-2\upsilon )(r_{,k}{n}_{,m}-r_{,m}{n}_{,k})\right\}$$

(10)

$$u_{km}^{\ast }(\mathit{x},\mathit{\xi })=\frac{1}{8\pi G(1-\nu )}[(3-4\upsilon )\ln \frac{1}{r}{\delta }_{km}+r_{,k}{r}_{,m}]$$

(11)

The integrals involving the kernel $t_{km}^{\ast }(\mathit{x},\mathit{\xi })$ are usually strongly singular, and in conventional IGABEM, it is usual to separate the jump term from the integral and produce a Cauchy principal-value integral [61], which is a rather complicated process.

In this paper, through the equation proposed by Liu [62,63], the jump term is transformed into the form of an integral, which has the following integral properties when the source point is at the boundary of the geometric model

$$C_{ij}(\mathit{\xi })=-{\displaystyle {\int }_{\Gamma }{t}_{ij}^{\ast }(\mathit{\xi },\mathit{x})d\Gamma (\mathit{x})}$$

(12)

Applying the above equation in Equation (9), the following weakly singular form of the BIE is obtained:

$${\displaystyle {\int }_{\Gamma }{t}_{ij}^{\ast }(\mathit{\xi },\mathit{x})[u_j(\mathit{x})-u_j(\mathit{\xi })]d\Gamma (\mathit{x})}={\displaystyle {\int }_{\Gamma }{u}^{\ast }(\mathit{\xi },\mathit{x})t_j(\mathit{x})d\Gamma (\mathit{x})}$$

(13)

in which, both of the integrals are weakly singular.

In the IGABEM, the shape functions, which are usually taken as general polynomials in the traditional BEM, are replaced with NURBS basis functions, thus maintaining all the geometric information of the original geometry model. And relationship $\mathit{x}=\mathit{x}(\xi )$ implies that in the parameter coordinates $\xi$, the displacements and tractions in Equation (13) are approximated by

$$u_j(\xi )={\displaystyle \sum _{l=1}^{p+1}{R}_{i,p}^{le}(\xi )d_j^{le}}$$

(14)

$$t_j(\xi )={\displaystyle \sum _{l=1}^{p+1}{R}_{i,p}^{le}(\xi )q_j^{le}}$$

(15)

where $d_j^{le}$ and $q_j^{le}$, denoting the displacement coefficients and traction coefficients of the lth control point of the eth element, respectively, are associated with a particular control point. Due to the localized support of the basis functions, only p + 1 basis functions are non-zero for each element.

In this paper, the boundary of the geometric model is divided into multiple independent NURBS curves with the endpoints of the C⁰ continuity. The boundary is discretized into a set of non-overlapping N_h patches

$$\Gamma =\underset{e=1}{\overset{N_h}{\cup }}{\Gamma }_h,{\Gamma }_i\cap {\Gamma }_j=0,i\ne j.$$

(16)

Each patch N_h is then discretized separately into a set of non-overlapping N_eh elements

$${\Gamma }_h=\underset{e=1}{\overset{N_{eh}}{\cup }}{\Gamma }_{eh},{\Gamma }_i\cap {\Gamma }_j=0,i\ne j.$$

(17)

The discretized boundary integral equation is then rewritten as

$$\begin{array}{l}{\displaystyle \sum _{h=1}^{N_h}{\displaystyle \sum _{e=1}^{N_e}{\displaystyle \sum _{l=1}^{p+1}\left[{\displaystyle {\int }_{-1}^{+1}{t}_{ij}^{\ast }\left({\mathit{\xi }}^{\prime },\mathit{\xi }(\widehat{\xi })\right)\left(R_{i,p}^{leh}\left(\widehat{\xi }\right)-R_{i,p}^{l\overline{e}h}\left({\widehat{\xi }}^{\prime }\right)\right)J^{eh}\left(\widehat{\xi }\right)d\widehat{\xi }}\right]}}} d_j^{leh}\\ ={\displaystyle \sum _{h=1}^{N_h}{\displaystyle \sum _{e=1}^{N_e}{\displaystyle \sum _{l=1}^{p+1}\left[{\displaystyle {\int }_{-1}^{+1}{u}_{ij}^{\ast }\left({\mathit{\xi }}^{\prime },\mathit{\xi }(\widehat{\xi })\right)R_{i,p}^{leh}\left(\widehat{\xi }\right)J^{eh}\left(\widehat{\xi }\right)d\widehat{\xi }}\right]}}} q_j^{leh}\end{array}$$

(18)

where $•{ }^{leh}$ is the quantity on the lth control point of the eth element of the hth patch, and $\overline{e}$ is the element in which the collocation point is located, and the Jacobian of transformation $J^{eh}\left(\widehat{\xi }\right)$ is given by [27]

$$J^{eh}(\widehat{\xi })=\frac{d{\Gamma }_{eh}}{d\xi }\frac{d\xi }{d\widehat{\xi }}=\sqrt{{\left(\frac{dx}{d\xi }\right)}^2+{\left(\frac{dy}{d\xi }\right)}^2}\frac{{\xi }_2-{\xi }_1}{2}$$

(19)

where ${\xi }_2$ and ${\xi }_1$ represent the upper and lower bounds of the element in the parameter space, respectively.

To obtain the linear system, a series of source points (collocation points) are selected on the boundary

$$\mathit{Hu}=\mathit{Gt}$$

(20)

where H is the coefficient matrix of the displacement vector u and G is the coefficient matrix of the traction vector t. Rearranging this set of equations by placing all unknown components on the left-hand side and known components on the right, the following can be written [27]:

$$\mathit{Ax}=\mathit{b}$$

(21)

where b is the vector containing the unknowns d and q. Mallardo and Ruocco [31] give a general method for the imposition of boundary conditions. In the present paper, the approach in the work by Simpson et al. [27,28] is used for the imposition of boundary conditions for simplicity.

3.2. Collocation Method

To obtain a system of equations in IGABEM, the collocation points must be given. In this paper, a new collocation method with no correction coefficient is proposed to move the points inside the patches. The collocation points are evaluated by the following formula

$${\xi }_a^{\prime }=\frac{{\xi }_a+{\xi }_{a+1}+\dots +{\xi }_{a+p+1}}{p+2} a=1,2,\dots ,n$$

(22)

The number of total collocation points is usually not equal to the number of total control points when this method is used because some control points may be shared by multiple patches. To ensure that the number of equations equals to the number of the control points, The method, proposed by Wang and Benson [32], for merging different equations related to the same control point to one equation is introduced

$$\frac{1}{2}({\overline{A}}_{\overline{m}1}{x}_1+{\overline{A}}_{\overline{m}2}{x}_2+\dots +{\overline{A}}_{\overline{m}n}{x}_n)={\overline{b}}_{\overline{m}}$$

(23)

where

$${\overline{A}}_{\overline{m}p}=A_{ip}+A_{\overline{j}p} p=1,2,\dots ,n$$

(24)

in which A_ip and A_jp represent the initial assembly matrix of two collocation points associated with corner points, which are directly added and combined.

3.3. Weakly Singular Numerical Integration Method

The coordinate transformation method proposed by Telles [60] is a simple and efficient method for the evaluation of weakly singular integrals. Compared with the element division method, this transformation technique uses fewer boundary elements and reduces the computational cost while maintaining higher accuracy. In the coordinate transformation method proposed by Telles, there exist two forms of transformation, i.e., second-order and third-order. The second-order transformation is suitable for dealing with the case where the singularities are located at the ends of the element. For the general case where the singularity is located inside the element, the element is split at the singular point and then the coordinate transformation method of the second order can be applied to both sides of the element separately. Conversely, the third-order form is able to handle singularities at arbitrary locations, but the corresponding equations are more complicated. Here, Telles’ method of order 2 is adopted to deal with the weakly singular integrals.

4. Optimization Algorithm

PSO Algorithm

PSO is an optimization technique proposed by Kennedy et al. [64] based on the study of bird flock foraging behavior. Individuals within a flock of birds search for the optimal destination while at the same time making the group find the collective optimal destination by continuously sharing information.

As shown in

Table 1. Algorithmic Principle Correspondence.

Flock of Birds	PSO
Bird	Particle
Forest	Solution space
Amount of food	Objective function value
Location of each bird	Particle position
The location with the most food	Global optimal solution

, the flock foraging behavior and the algorithm principle are corresponded. In the PSO algorithm, the initial value is a set of particles randomly distributed within the D-dimensional solution space, each particle corresponds to an objective function value, and the position of the ith particle at the kth iteration is:

$$x_i^k=(x_{i1}^k,x_{i2}^k,\dots ,x_{iD}^k)$$

(25)

The position update vector of a particle is called velocity, and each particle has its own position and velocity. The individual historical optimum in an iterative step is called the individual extreme value (x_pbest), and the group optimum in this step is called the group extreme value (x_gbest). In the PSO algorithm, each particle is regarded as a potential solution in the problem search space, and these particles move iteratively in the solution space and adjust their flight paths by tracking their own individual optimal positions as well as the global optimal position of the whole group. The particle position and velocity iteration formulas are [65]

$$v_i^{k+1}=w{v}_i^k+c_1{r}_1^k(x_{pbest}-x_i^k)+c_2{r}_2^k(x_{gbest}-x_i^k)$$

(26)

$$x_i^{k+1}=x_i^k+v_i^{k+1}$$

(27)

where k is the number of iterations; w is the inertia weight, which represents the dependence of the particle on its previous motion state; c₁ and c₂ are learning factors, which, respectively, adjust the maximum step towards the individual best particle and the global best particle; and $r_1^k$ and $r_2^k$ are two random numbers between 0 and 1.

A fitness function is a function, incorporating the objective function and constraints to assess the quality of solutions, defined as follows [59]

$$fitness=f\left(x\right)+\zeta \times \max \left(g_j,0\right) +{\zeta }^{\prime }\times \left|h_i\right|$$

(28)

where f(x) is the objective function, g_j and h_i, respectively, represent the inequality constraints and equality constraints, and $\zeta$ and ${\zeta }^{\prime }$ are the penalty factors for the constraints. In this study, the equality constraints are not considered, i.e., h_i = 0. The fitness functions are given in every optimization design example.

5. Numerical Examples of Multi-Patch WSIGABEM

In order to verify the accuracy of the algorithm, the multi-patch WSIGABEM is compared with analytical solutions and the traditional IGABEM. Starting from the simple case of a fillet plate model subjected to uniaxial tension, the algorithm is further applied to a more complex plate with holes. In order to investigate the advantages and disadvantages of the collocation methods in the multi-patch WSIGABEM, the newly proposed collocation point method, i.e., the no-coefficient collocation method, is compared with the improved collocation method with coefficient by Wang et al. [32], which has optimal coefficients for different n (number of mesh refinements), d (number of Gauss points), and p (number of basis function orders).

The following definition is used to calculate the relative L₂ error norm in displacements on the boundary [27]:

$$e_{L_2}=\frac{\left|\right|\mathit{u}-{\mathit{u}}_{\mathit{e}\mathit{x}}|{|}_{L_2}}{\left|\right|{\mathit{u}}_{\mathit{e}\mathit{x}}|{|}_{L_2}}$$

(29)

5.1. A Plate under Uniaxial Tension

A plate under uniaxial tension is considered. Due to the symmetry of the problem, a quarter of the plate, shown in Figure 1, is sufficient to represent the mechanical behavior of the entire plate. When performing the analysis, the corresponding symmetric boundary conditions are applied. This approach is common in engineering analysis and can greatly improve the efficiency [1]. Figure 2 shows the different meshes used in the following numerical calculations, in which the order p is 2 or 3.

Figure 3 shows the numerical errors of the displacement on the boundary calculated by the multi-patch WSIGABEM with the collocation points suggested by Wang et al. [32]. As shown in Figure 3a, when p = 2 and d = 4, the optimal value of β is 0.6 for the different meshes. As shown in Figure 3b, when p = 2 and d = 8, the optimal value of β is 0.3 for the different meshes. When p = 3 and d = 4, as shown in Figure 3c, the optimal value of β is 0.2; and when p = 3 and d = 8, as shown in Figure 3d, the optimal value for β is 0.45. It seems that the meshes have no obvious influence on the value of the optimal coefficient β. The optimal choice of the coefficient β is empirical. It is influenced by the number of Gauss points, the order of the curve, and the refinement of the mesh. From Figure 3a,b, or Figure 3c,d, it is obvious that the optimal value of β is different when the number of the Gauss points is 4 and 8. A comparison of Figure 3a with Figure 3c, or Figure 3b with Figure 3d, shows that the optimal value of β is different when the order of the curves is 2 and 3. All four subfigures in Figure 3 show the insensitivity of the coefficient β toward the mesh refinement.

Figure 4 is a comparison of the results by the IGABEM [27] and the multi-patch WSIGABEM with the collocation points suggested by Wang et al. [32]. The multi-patch WSIGABEM has a higher accuracy when the number of the Gauss points is low and the degree of freedom is small when the value of coefficient β is optimal.

Figure 5 shows the numerical errors of the results calculated by the multi-patch WSIGABEM with three different methods for the generation of the collocation points. The three methods are the no-coefficient collocation point method and the improved collocation point method with coefficient β set to 0.3 and 0.6. As shown in Figure 5a, with the same refinement, n = 4, and the same order of basis functions, p = 2, the no-coefficient collocation point method has high accuracy and good convergence, and its error is between the other two methods. In addition, the relationship between the error and the degree of freedom is shown in Figure 5b, which shows that the errors of the three methods are in the same order. Therefore, the new collocation method has quite similar accuracy to the improved collocation method by Wang et al. [32], while it locates the collocation points inside the patch without any correction.

5.2. An Infinite Plate with a Hole

The problem of a hole within an infinite plate subjected to unidirectional stretching has a known analytical solution [66]. Given the symmetry of the structure and boundary conditions, the problem can be reduced to analyze a quarter plate. To facilitate the numerical solution, the infinite-sized plate problem is equated to a finite-sized plate. A finite rectangular plate is intercepted while replacing the boundary force with an analytical traction at the boundary. The model of the problem is given in Figure 6, in which one-way displacement constraints are imposed on the lower and right edges, while analytic tractions are imposed on the upper and left edges. When the NURBS curves are quadratic and cubic (p = 2 and p = 3), the meshes of the boundary are shown in Figure 7.

Figure 8 shows the numerical results calculated by the multi-patch WSIGABEM for the refinement, the curve order, and the number of Gaussian points. The optimal values of the coefficient β are shown in Figure 8 and listed in

Table 2. The optimal correction coefficients for different d and n, when p = 2.

p = 2	mesh3	mesh4	mesh5
d = 4	0.15	0.15	0.15
d = 8	0.1	0.1	0.1

and

Table 3. The optimal correction coefficients for different d and n, when p = 3.

p = 3	mesh3	mesh4	mesh5
d = 4	0.75	0.8	0.8
d = 8	0.2	0.55	0.6

. By comparing the data in the tables, it is obvious that the number of Gauss points and the order of the basis function are the main influencing factors. As shown in Figure 8c,d, The relationship between the correction coefficient and different mesh refinements is different. The optimal value of β changes with the number of mesh refinements, which is different from the conclusion of the square plate case.

This difference may be related to the complexity of the problem. The optimal values of the correction coefficients β are also inconsistent in different cases. As shown in Figure 3 and Figure 8, different values of correction coefficients may lead to an order-of-magnitude difference in computational accuracy. Therefore, in practical application, it is necessary to select the appropriate coefficient for the specific analysis to achieve the best calculation accuracy.

Figure 9 demonstrates that the multi-patch WSIGABEM has a higher accuracy when the number of the Gauss points is low and the degree of freedom is small when the value of coefficient β is optimal.

Figure 10 shows a comparison with different collocation point methods in the multi-patch WSIGABEM. According to Table 2 and Table 3, the coefficients β are chosen to be 0.1 and 0.15. Although the improved collocation point method seems to be more accurate, it depends on the choice of the coefficient. The newly proposed collocation point method has comparable accuracy, while it has the advantage in the implementation.

6. Optimization Design Based on the Multi-Patch WSIGABEM

In this section, three optimal designs of the shapes of the fillet corner, spanner, and arch bridge based on the multi-patch WSIGABEM and PSO are presented for the verification of the accuracy and effectiveness of this approach. For simplicity, the materials in the examples are assumed to be the same, with Young’s modulus E = 10⁷ units and Poison’s ratio v = 0.3, while the yield strengths are different. The number of Gaussian points is 4. For better accuracy, the order of the curves is p = 3, and the mesh is refined to be mesh7 following the same procedure as in Figure 2 and Figure 7. The von Mises stresses on the boundaries are calculated by the stress recovery method explained in Appendix A. The learning factors c₁ and c₂ are 1.5, the maximum flight speed v is 1/5 of the search space, and the inertia weights in the PSO method are improved by:

$$w=w_{\max }-(w_{\max }-w_{\min })\frac{i}{i_{\max }}$$

(30)

where the maximum weight $w_{\max }$ and the minimum weight $w_{\min }$ are, respectively, set to be 0.6 and 0.2 in the three examples, i is the step of the current iteration, and i_max is the maximum number of iterations. The weights $w$ are adjusted by using the adaptive linear method, which enables the algorithm to have a strong global search ability in the early stage and a better local convergence speed in the later stage. These parameters are chosen to ensure the convergence and efficiency of the optimization algorithm. The accuracy of the algorithm is verified by the comparison of the optimal shape of the fillet corner by the present method and by the existing method [57,59]. Then, the shape optimizations of the spanner and the arch bridge are investigated. The coordinates and weights of the control points of the shapes before the optimizations are given in Appendix B.

6.1. The Shape Optimization of the Fillet Corner

As shown in Figure 11a, a normal traction t = 100 units is prescribed at the blue line, resulting in a stress concentration at point N. The shape optimization of the fillet corner is conducted to find the minimum area while the von Mises stress is no greater than the yield strength of 125 units. The symmetric condition is prescribed on the left and bottom sides, denoted as red lines since the shape in Figure 11 is only one-quarter of the whole fillet. The design parameters are chosen to be the coordinates and weights of the control points A and B, denoted as red points in the figure. The shape of the corner changes when the design parameters change. During the optimization, point A is assumed to be at the top left of point B.

The initial values of the design parameters are randomly distributed in the ranges listed in

Table 4. The ranges and the values of the design parameters before and after the optimization for the shape optimization of the fillet.

Design Parameters	Lower Bound	Upper Bound	Before Optimization	After Optimization
A	(9, 4.5)	(15, 9)	(11.17, 1.5)	(9.00, 8.00)
B	(9, 4.5)	(15, 9)	(13.33, 6)	(9.00, 4.50)
wA	0.1	10	1	2.49
wB	0.1	10	1	1.74

. The fitness function is defined as follows:

$$fitness=f\left(x\right)+\zeta \times \max \left(vonMise{s}_{\max }-125,0\right)$$

(31)

where $f\left(x\right)$ is the objective function, representing the area of the part, $vonMise{s}_{\max }$ represents the maximum von Mises stress value on the boundary, and $\zeta$ is the penalty factor, which is chosen to be 1 after some trials. At the initial stage, there are 20 sets of design parameters and the maximum number of iterations is 20. After the optimization procedure, the shape of the fillet is as shown in Figure 11b. The area of the fillet with the minimum value of the fitness function on every iteration step is as shown in Figure 12. The last column of Table 4 gives the values of the design parameters after the shape optimization. To speed up the iteration process, one set of the design parameters at the initial stage is chosen to be the same as that shown in Figure 11a, i.e., A(11.17, 7.5), B(13.33, 6), and weights w_A = w_B = 1.

Figure 13a,b show von Mises stresses on the boundary before and after the optimization, respectively, with the first and last knot values of the closed curves corresponding to the coordinates starting from the point (0, 0) counterclockwise to the end of the point (0, 0). A similar shape optimization was conducted by Lian et al. [57] and Sun et al. [59], in which only the vertical coordinates of the control points served as the design parameters. The optimal area predicted by the former works was about 138 units, while in the present design, all the coordinates and weights of the control points serve as the design parameters. The optimal area by the present method is 132.46 units, which is less than 138 units. Due to the PSO algorithm, the complicated analysis of sensitivity is not necessary, making the realization process simpler and easier.

6.2. The Shape Optimization of the Spanner

Figure 14 shows the optimal design of the shape of the spanner. The boundaries denoted by the red lines are fixed and the shear traction t = 10 units is prescribed at the blue line. The vertical coordinates and the weights of the red points are the design parameters to be optimized. Since the control points are symmetrically distributed, only the top four control points are set as the design parameters. In addition, for simplicity, the weights of points A, D remain the same during the optimization. In summary, there are six design parameters, which are the vertical coordinates y₁, y₂, y₃, and y₄ of points A, B, C, and D, and the weights w₁ and w₂ of points B and C. To speed up the iteration process, one set of the design parameters at the initial stage is chosen to be A(−2, 1.5), B(−1, 3), C(1, 1), B(5, 1), and weights w₁ = w₂ = 2.

The shape optimization of the spanner is conducted to find the minimum area while the von Mises stress is no greater than the yield strength of 825 units, and the fitness function is defined as follows

$$fitness=f\left(x\right)+\zeta \times \max \left(vonMise{s}_{\max }-825,0\right)$$

(32)

where $f\left(x\right)$ is the area of the spanner and the penalty factor $\zeta =10$. At the initial stage, there are 20 sets of design parameters. The maximum number of iterations is set to be 20. After the optimization procedure, the shape of the spanner is as shown in Figure 14b, while the trajectory of the area in the iteration procedure is as shown in Figure 15. The range of design parameters and the optimized design parameter values are shown in

Table 5. The ranges and the values of the design parameters before and after the optimization for the shape optimization of the spanner.

Design Variable	Lower Bound	Upper Bound	Before Optimization	After Optimization
A	(−2, 1.1)	(−2, 3)	(−2, 3)	(−2, 1.18)
B	(−1, 1)	(−1, 3)	(−1, 3)	(−1, 2.33)
C	(1, 0)	(1, 3)	(1, 3)	(1, 1.17)
D	(5, 1)	(5, 3)	(5, 3)	(5, 1.00)
wB	0.1	10	1	2.93
wC	0.1	10	1	1.52

. The area of the spanner is reduced from 68.43 units to 25.96 units, which is only 37.93% of the original shape.

6.3. The Shape Optimization of the Arch Bridge

The shape optimization of the arch bridge is shown in Figure 16. The fixed constraint is applied at the red line and the normal traction t = 10 units is prescribed at the blue line, which is at the center of the bridge. The design parameters are related to the red points. Due to the symmetry, only points A and B in the figure are optimized. Since the weight of point B is the main factor affecting the shape of the boundary, the weight of point A remains the same during the optimization. The design parameters are the horizontal coordinates of point A and the coordinates and weight of point B. To speed up the iteration process, one set of the design parameters at the initial stage is chosen to be A(4, 0), B(4, 1), and weight w_B = $1/\sqrt{2}$.

The shape optimization of the Arch bridge is conducted to find the minimum area while the von Mises stress is no greater than the yield strength of 50 units, and the fitness function is defined as follows

$$fitness=f\left(x\right)+\zeta \times \max \left(vonMise{s}_{\max }-50,0\right)$$

(33)

In this example, the penalty factor $\zeta =1$, the number of particles is taken as 20, and the maximum number of iterations is 20. After the optimization procedure, the shape of the arch bridge is as shown in Figure 16b, while the trajectory of the area iteration process is as shown in Figure 17. The range of design parameters and the values of the design parameters before and after the optimization are shown in

Table 6. The ranges and the values of the design parameters before and after the optimization for the shape optimization of the arch bridge.

Design Variable	Lower Bound	Upper Bound	Before Optimization	After Optimization
A	(0, 0)	($4\sqrt{3}$, 0)	($2\sqrt{3}$, 0)	(3.20, 0)
B	(0, 0)	($2\sqrt{3}$, 2)	($\sqrt{3}$, 0)	(1.81, 1.34)
wB	0.1	10	1	1.85

. The area is reduced from 20.78 units to 14.73 units, which is only 70.88% of the original.

7. Conclusions

A multi-patch WSIGABEM is implemented in the present paper. A new formula is proposed for the generation of the collocation points to keep the points off the discontinuity boundaries. The utility of the multi-patch WSIGABEM is verified by comparison with analytical solutions and the traditional IGABEM. The multi-patch WSIGABEM with the newly proposed formula for the generation of the collocation points has a similar accuracy to that of the improved Greville abscissae method. However, the improved Greville abscissae method requires experience in the choice of coefficients, while the newly proposed formula is more straightforward. The new collocation method is suggested to simplify the implementation of the multi-patch WSIGABEM.

In addition, the combination of the multi-patch WSIGABEM and the PSO algorithm yields a novel shape optimization algorithm, in which the coordinates and weights of the control points are used as the design parameters. As proved by three examples, the method effectively improves structural performance, reduces weight, and lowers the cost, and thus provides a feasible algorithm in the field of structural shape optimization.