A Gradient-Free Topology Optimization Strategy for Continuum Structures with Design-Dependent Boundary Loads

: In this paper, the topology optimization of continuum structures with design-dependent loads is studied with a gradient-free topology optimization method in combination with adaptive body-ﬁtted ﬁnite element mesh. The material-ﬁeld series-expansion (MFSE) model represents the structural topology using a bounded material ﬁeld with speciﬁed spatial correlation and provides a crisp structural boundary description. This feature makes it convenient to identify the loading surface for the application of the design-dependent boundary loads and to generate a body-ﬁtted mesh for structural analysis. Using the dimension reduction technique, the number of design variables is signiﬁcantly decreased, which enables the use of an efﬁcient Kriging-based algorithm to solve the topology optimization problem. The effectiveness of the proposed method is demonstrated using several numerical examples, among which a design problem with geometry and contact nonlinearity is included.


Introduction
During the past three decades, topology optimization has been widely applied to determine the optimal material distribution of various structural and multidisciplinary design problems [1][2][3][4][5][6][7]. Compared with structural topology optimization with invariant loads, the main feature of the optimization problems with design-dependent loads is that the load boundaries keep evolving during the optimization procedures. According to the load direction, design-dependent loads can be generally classified into two categories: (i) the load direction is always perpendicular to the loading surface, such as the fluid pressure load on structures in hydrostatics problems; (ii) the load direction is fixed, such as the soil pressure load on civil structures. The main issues of topology optimization with design-dependent loads lie in (i) identifying the surface on which a load is applied during an optimization process and (ii) the difficulty in obtaining the accurate sensitivity information of design-dependent boundary loads.
In terms of identifying the evolving loading surface in topology optimization, the commonly used approaches are mainly based on the element-density method. Typically, Hammer and Olhoff [8] proposed a parameterized isovolumetric density curve to represent the loading surface, where a load is distributed to nodes using an interpolation function for FE analysis. To avoid the appearance of an invalid loading surface, Du and Olhoff [9] further suggested a modified isoline technique, where the loading surface is determined based on the successive isoline information. Fuchs and Shemesh [10] introduced a parametric loading surface defined by the polynomial curve or Bezier curve, and the surface parameters were used as additional design variables of the optimization problem. Zhang et al. [11] proposed an element-based search scheme, and they used it to identify loading surfaces according to The original design is shown in Figure 1a, and a pressure force P is applied on the top of the original structure. As the optimization iterations continue, the original structure changes, and the locations and directions of the applied pressure load P changes accordingly (see Figure 1b). In this paper, the topology optimization aimed at minimizing the compliance of a structure subjected to design-dependent pressure loads with given material usage, and The original design is shown in Figure 1a, and a pressure force P is applied on the top of the original structure. As the optimization iterations continue, the original structure changes, and the locations and directions of the applied pressure load P changes accordingly (see Figure 1b).
In this paper, the topology optimization aimed at minimizing the compliance of a structure subjected to design-dependent pressure loads with given material usage, and the formulation is as follows: find χ min C = P(χ) T u(χ) s.t. K(χ)u = P(χ) where C is the structure mean compliance and the 0/1 discrete function χ denotes the structural topology, i.e., χ = 0 signifies a void material, while χ = 1 signifies a solid material. K is the stiffness matrix. u is the structural nodal displacement vector obtained from the linear or nonlinear equilibrium equation K(χ)u = P(χ). V * is the constant volume of the material constraint.

Bounded Material Field Definition
Here, a bounded material-field function ϕ(x) ∈ [−1, 1] with a certain degree of correlation behavior was defined in the design domain x ∈ Ω x to describe the structural topology χ. The mathematical description of the structural topology with the material-field function was represented as follows: where d z is a given threshold which satisfies −1 < d z ≤ 0.
For the material-field function ϕ(x), the spatial correlation at the two points closely depends on their spatial distance, and an exponential model was used to describe the correlation as follows: where x 1 and x 2 are any two points in the design domain x ∈ Ω x , and L c is the correlation length of the material field. According to Equation (3), the correlation length L c is an important factor for determining the correlation between the two points. Therefore, the correlation length L c determined the spatial fluctuation degree of the material field. When the structural topology was represented by the material field ϕ(x), the value of the correlation length L c could be used to control the degree of detail of the structures [24]. For a 2D rectangular design domain with dimensions l a and l b , the correlation length was set to L c = 0.3 · min(l a , l b ) in this study.

Reduced Series Expansion of the Material Field
Based on the theory of bounded fields [37][38][39], the design domain was discretized into a finite number of regularly arranged material-field observation points x i (i = 1, 2, . . . , N P ), and the series expansion of the material-field ϕ(x) was expressed by x N P is a vector, that contains the correlation between x and all the material-field points, and ξ = ξ 1 ξ 2 . . . ξ N P T is a vector of uncorrelated MFSE coefficients. α j and ψ j are the eigenvalues and eigenvectors of the correction matrix R, which is expressed as follows: The eigenvalues α j of R are in descending order. The greater the eigenvalues, the greater the contribution to the corresponding expansion item in Equation (4). Therefore, considering the question of improving the calculation efficiency without losing sufficient accuracy, it is efficient to truncate Equation (4) to only a few main M items (M << N P ). The truncation criterion is as follows: where ε is a small number. In this paper, we set ε = 0.001. After the truncation of Equation (4), the material-field ϕ(x) could be approximately expressed in the vector form, which is where ξ = {ξ 1 ξ 2 . . . ξ M } T contains the M uncertain coefficients. Λ is the diagonal matrix that consists of M eigenvalues α j , and Ψ consists of M eigenvectors ψ j . According to the range of the material-field ϕ(x, ξ) ∈ [−1, 1] and Equation (7), we could obtain By squaring the Equation (8), we get By introducing the matrix

Identification of the Loading Surface
For topology optimization problems with design-dependent loads, identifying the loading surface is one of the major tasks. In the MFSE method, it is very convenient to identify the loading surface as where the non-load-carrying parts represent the void material, and the load-carrying parts represent the solid material. In the boundary and topology description of Equation (11), the different realizations of the bounded material field ϕ(x, ξ) and different thresholds d z represent different loading surfaces. In this study, the value of the threshold d z is not a fixed value and the maximum value of d z is 0. The minimum value of d z depends on the volume constraint V * in the optimization formula Equation (1). Generally speaking, the smaller the volume constraint V * , the smaller the minimum value of d z . The value of d z gradually increased with the optimization process.
The schematic diagram for the loading surface and the body-fitted mesh based on the MFSE is shown in Figure 2. As shown in Figure 2a, we assumed that the design domain is rectangular with a dimension of l a × l b . Initial pressure forces were applied on the top edge of the design domain, as shown in Figure 2b and the final loading surface is illustrated in Figure 2c. The whole process mainly consists of three steps as follows:  and that the position of the material-field point is represented by 1 2 , , ..., According to the correlation formulation in Equation (3), the correction ma- Step 1: Represent the Material-Field Function Using Series Expansion Suppose that there are N P observation points evenly distributed in the rectangular design domain x ∈ Ω x and that the position of the material-field point is represented by x 1 , x 2 , . . . , x N P . According to the correlation formulation in Equation (3), the correction matrix R containing the correlation information of all the material-field points can be determined. When a set of uncorrelated MFSE coefficients (design variables) ξ * = {ξ 1 * ξ 2 * . . . ξ M * } T is given, the material-field function value ϕ(x, ξ * ) in the design domain could be determined using Equation (7). In Figure 2a, a material-field function ϕ(x, ξ * ) was given, where its upper and lower boundaries are between −1 and 1. The green plate is the given threshold d z , and it was used to cut the material-field function ϕ(x, ξ * ) and divide the design domain into load-carrying parts and non-load-carrying parts according to Equation (11).
Step 2: Determine the Load-Carrying Parts and Update the Loading Surface After cutting the material-field function ϕ(x, ξ * ) using the given threshold d z , the contours of the function value ϕ(x, ξ * ) in the design domain are obtained. For the case ϕ(x, ξ * ) ≥ d z , the load-carrying parts are denoted by the red part. The white part in Figure 2c represents the non-load-carrying part that satisfies ϕ(x, ξ * ) < d z . There is no doubt that the loads were applied to the load-carrying parts. The detailed process of updating the loading surfaces of any structure: (i) When the initial pressure is initially applied on the front edge of the design domain (e.g., Figure 2b), we choose N discrete points X m (m = 1, 2, . . . , N) uniformly distributed on the initial loading surface as shown in Figure 3a, and we define the distance between adjacent discrete points X m and X m+1 is d. (ii) Judging the material-field value from each point X m along the pressure direction and searching out the positionx m that first satisfies ϕ(x m , ξ * ) ≥ d z (represented by the yellow dots in Figure 3a). It should be pointed out that if ϕ(x m , ξ * ) ≥ d z does not exist in the design domain, the loading pointx m is marked in the back edge of the design domain.
(iii) Connecting all the pointsx m to spline curves which are obtained by the cubic spline interpolation. Based on the above process, the blue curves in Figure 3 for two typical cases are identified as the loading surfaces.   During the proposed gradient-free optimization process, the sample is consid invalid if the loading surface is discontinuous as shown in Figure 3b. These invalid s ples are easily identified by judging whether the distance between two jacent discrete points ˆ x and ˆ x on the loading surface is greater than a given v During the proposed gradient-free optimization process, the sample is considered invalid if the loading surface is discontinuous as shown in Figure 3b. These invalid samples are easily identified by judging whether the distance dis(x m ,x m+1 ) between two adjacent discrete pointsx m andx m+1 on the loading surface is greater than a given value D= 5 · d.
Step 3: Apply the Pressure Loads and Re-Meshing for the Finite Element Analysis According to the contour obtained in Step 2, the design domain was cut into several parts: load-carrying parts and non-load-carrying parts. To avoid mesh generation failure for the cutting design domain, it was reasonable to mesh with irregular triangle elements, and the process of meshing is shown in Figure 4. The pressure loads were then applied to the nodes of the triangular elements on the loading surface, as shown in Figure 2d. Suppose that there were e N triangle elements (e 1 , . . . , e N ) in the design domain and that the barycenter of the triangle elements is represented by (x e 1 , . . . , x e N ), as shown in Figure 2d by the black points. According to the mapping of the material-field function value ϕ(x, ξ * ), we could determine the material-field function value ϕ(x e 1 , ξ * ), . . . , ϕ(x e N , ξ * ) at the elemental barycenter x e 1 , . . . , x e N . Then, based on the material interpolation function, the material property of each element in Figure 2d was determined to complete the finite element analysis.
Symmetry 2021, 13, 1976 8 of 19 for the cutting design domain, it was reasonable to mesh with irregular triangle elements, and the process of meshing is shown in Figure 4. The pressure loads were then applied to the nodes of the triangular elements on the loading surface, as shown in Figure

Topology Optimization Formulation for the Structures with Design-Dependent Boundary Loads
According to Equation (2) and the material-field function ( ) , it is reasonable to use the SIMP interpolation function to map the material-field function ( ) , ϕ x ξ to the elastic modulus of the structure. Here, we assumed that there were N e finite elements, and the interpolation function of the elemental elastic modulus was shown as follows: where 0 E is the elastic modulus of fully solid materials, min E is a small positive value to avoid singularity, and it was assumed that is the Heaviside function projection [40] of the material-field function ( ) , ϕ x ξ , which was used to obtain a clear topology structure, where the function is

Modeling the initial design domain
Obtain the material field function according to Eq. (7) Output the material-field contours by given the threshold

Topology Optimization Formulation for the Structures with Design-Dependent Boundary Loads
According to Equation (2) and the material-field function ϕ(x, ξ) ∈ [−1, 1], it is reasonable to use the SIMP interpolation function to map the material-field function ϕ(x, ξ) to the elastic modulus of the structure. Here, we assumed that there were e N finite elements, and the interpolation function of the elemental elastic modulus was shown as follows: where E 0 is the elastic modulus of fully solid materials, E min is a small positive value to avoid singularity, and it was assumed that E min = 10 −6 E 0 . ϕ(x, ξ) is the Heaviside function projection [40] of the material-field function ϕ(x, ξ), which was used to obtain a clear topology structure, where the function is where β denotes the smoothing parameter. According to Equation (7), the minimum compliance topology optimization problem using the MFSE method was expressed as follows: where V e denotes the eth elemental volume, and the design variables are the MFSE coefficients ξ = {ξ 1 , ξ 2 , . . . , ξ M } T in the ξ-space. Note that the elemental relative density could be determined using ρ e (ξ) = 1+ ϕ(x e ,ξ) 2 .

Sequential Kriging-Based Optimization Algorithm
In this study, we used the sequential optimization algorithm [24] based on the Kriging surrogate model to complete the MFSE-based topology optimization problem (15) with design-dependent loads. According to the algorithm, in each sub-domain Ω k , the topology optimization model in Equation (14) needs to be converted to the unconstrained form, which is: where p 0 is the penalty multiplier.
The sequential Kriging-based optimization algorithm consists of several suboptimization problems, where each sub-optimization problem contains two steps: (i) choosing initial samples using Latin hypercube sampling to construct a surrogate model [41]; (ii) adding post samples using the EI criterion to obtain a smaller objective function value. The number of the post sampling N s is 30 in this paper.
To improve the accuracy of the surrogate model and the searching efficiency as well as reduce the number of invalid loading surfaces during sampling, it is necessary to use the self-adaptive strategy of adjusting the design domain to deal with the topology optimization problem with design-dependent pressure loads in the Kriging-based algorithm. For the specific implementation processes of the self-adaptive strategy, see literature [24]. In this paper, the structural samples of the invalid loading surface will produce poor performance, and it should be pointed out that the Kriging-based optimization algorithm can effectively avoid the omission of the invalid loading surface (judged by dis(x m ,x m+1 ) > D) in Section 3.3.
The flow chart of gradient-free topology optimization with design-dependent boundary loads is illustrated in Figure 5. The convergence criterion of the optimization procedure is that the change of the objective function in two consecutive iterations is less than 0.01.

Numerical Examples
To validate the effectiveness of the proposed method for the topology optimization with design-dependent pressure loads, we here presented five examples, including an example with contact supports with friction. In each example, Young's modulus and Poisson's ratio of the material were 0 200 GPa E = , and , respectively. The initial value of the smoothing parameter β was zero, and it was gradually increased to 20 during the optimization process. The plane stress assumption was considered for all the examples. The finite element re-meshing and analysis were completed in the commercial package ABAQUS. The number of design variables in Equation (14) is 38 for examples 1, 2, and 3, 53 for example 4, and 63 for example 5.

Example 1
In this example, a bridge-like rectangular structure with sides 120 mm x l = and 60 mm y l = was considered, as shown in Figure 6. The structure was clamped at the two bottom corners and subjected to an initial pressure =10 N/mm P from its top edge. The directions of the pressure loads were fixed, and the allowable material volume fraction was 50%.

Numerical Examples
To validate the effectiveness of the proposed method for the topology optimization with design-dependent pressure loads, we here presented five examples, including an example with contact supports with friction. In each example, Young's modulus and Poisson's ratio of the material were E 0 = 200 GPa, and υ = 0.3, respectively. The initial value of the smoothing parameter β was zero, and it was gradually increased to 20 during the optimization process. The plane stress assumption was considered for all the examples. The finite element re-meshing and analysis were completed in the commercial package ABAQUS. The number of design variables in Equation (14) is 38 for examples 1, 2, and 3, 53 for example 4, and 63 for example 5.

Example 1
In this example, a bridge-like rectangular structure with sides l x = 120 mm and l y = 60 mm was considered, as shown in Figure 6. The structure was clamped at the two bottom corners and subjected to an initial pressure P = 10 N/mm from its top edge. The directions of the pressure loads were fixed, and the allowable material volume fraction was 50%. The optimized design is shown in Figure 7 and its compliance is In Figure 7, the red part is the load-carrying part, where the new loading surface remained on the top edge of the load-carrying part. After cutting the design domain with the contour, we meshed the design domain with the triangle element shown in Figure 8. Based on the interpolation function in Equation (11), the elemental relative density in Figure 8 could be obtained. Figure 9 shows the optimized solution for the bridge-like structure example obtained with a gradient-based method by Du [9]. Compare the optimized topology in Figures 7 and 9, it can be seen that these two results are very similar. The two optimized solutions consist of a global single arch without small columns, and the curvature center of the loading surface is located above the top surface of the structure. In addition, for the optimal solution of reference [9] (Figure 9), the same load and boundary conditions as in Example 1 are applied to the structure, and the compliance of the structure is 118.82 Nmm C = , which is almost the same as the compliance of the optimal structure in this paper (    The optimized design is shown in Figure 7 and its compliance is C = 118.29 Nmm. In Figure 7, the red part is the load-carrying part, where the new loading surface remained on the top edge of the load-carrying part. After cutting the design domain with the contour, we meshed the design domain with the triangle element shown in Figure 8. Based on the interpolation function in Equation (11), the elemental relative density in Figure 8 could be obtained. Figure 9 shows the optimized solution for the bridge-like structure example obtained with a gradient-based method by Du [9]. Compare the optimized topology in Figures 7 and 9, it can be seen that these two results are very similar. The two optimized solutions consist of a global single arch without small columns, and the curvature center of the loading surface is located above the top surface of the structure. In addition, for the optimal solution of reference [9] (Figure 9), the same load and boundary conditions as in Example 1 are applied to the structure, and the compliance of the structure is C = 118.82 Nmm, which is almost the same as the compliance of the optimal structure in this paper (C = 118.82 Nmm vs. C = 118.29 Nmm). The optimized design is shown in Figure 7 and its compliance is In Figure 7, the red part is the load-carrying part, where the new loading surface remained on the top edge of the load-carrying part. After cutting the design domain with the contour, we meshed the design domain with the triangle element shown in Figure 8. Based on the interpolation function in Equation (11), the elemental relative density in Figure 8 could be obtained. Figure 9 shows the optimized solution for the bridge-like structure example obtained with a gradient-based method by Du [9]. Compare the optimized topology in Figures 7 and 9, it can be seen that these two results are very similar. The two optimized solutions consist of a global single arch without small columns, and the curvature center of the loading surface is located above the top surface of the structure. In addition, for the optimal solution of reference [9] (Figure 9), the same load and boundary conditions as in Example 1 are applied to the structure, and the compliance of the structure is 118.82 Nmm C = , which is almost the same as the compliance of the optimal structure in this paper (    The optimized design is shown in Figure 7 and its compliance is In Figure 7, the red part is the load-carrying part, where the new loading surface remained on the top edge of the load-carrying part. After cutting the design domain with the contour, we meshed the design domain with the triangle element shown in Figure 8. Based on the interpolation function in Equation (11), the elemental relative density in Figure 8 could be obtained. Figure 9 shows the optimized solution for the bridge-like structure example obtained with a gradient-based method by Du [9]. Compare the optimized topology in Figures 7 and 9, it can be seen that these two results are very similar. The two optimized solutions consist of a global single arch without small columns, and the curvature center of the loading surface is located above the top surface of the structure. In addition, for the optimal solution of reference [9] (Figure 9), the same load and boundary conditions as in Example 1 are applied to the structure, and the compliance of the structure is

Nmm C =
, which is almost the same as the compliance of the optimal structure in this paper (      The iteration history of the proposed MFSE-based gradient-free method with adaptive body-fitted meshes for solving the topology optimization problem with design-dependent loads was plotted in Figure 10. The objective function converged to a stable value during the iteration process while adjusting the design domain. The reason for the appearance of the fluctuant objective function is due to the Latin hypercube sampling in each sub-optimization.

Example 2
In this example, the optimization of a cover-like structure with the directions of the pressure loads perpendicular to the loading surface was considered. Figure 11    The iteration history of the proposed MFSE-based gradient-free method with adaptive body-fitted meshes for solving the topology optimization problem with design-dependent loads was plotted in Figure 10. The objective function converged to a stable value during the iteration process while adjusting the design domain. The reason for the appearance of the fluctuant objective function is due to the Latin hypercube sampling in each sub-optimization.  The iteration history of the proposed MFSE-based gradient-free method with adaptive body-fitted meshes for solving the topology optimization problem with design-dependent loads was plotted in Figure 10. The objective function converged to a stable value during the iteration process while adjusting the design domain. The reason for the appearance of the fluctuant objective function is due to the Latin hypercube sampling in each sub-optimization.

Example 2
In this example, the optimization of a cover-like structure with the directions of the pressure loads perpendicular to the loading surface was considered. Figure 11

Example 2
In this example, the optimization of a cover-like structure with the directions of the pressure loads perpendicular to the loading surface was considered. Figure 11 shows the design domain, boundary conditions, and initial pressure loads P = 10 N/mm. The sides of the rectangular design domain are l x = 120 mm and l y = 60 mm, and the allowable material volume fraction is 30%.  The iteration history of the proposed MFSE-based gradient-free method with adaptive body-fitted meshes for solving the topology optimization problem with design-dependent loads was plotted in Figure 10. The objective function converged to a stable value during the iteration process while adjusting the design domain. The reason for the appearance of the fluctuant objective function is due to the Latin hypercube sampling in each sub-optimization.

Example 2
In this example, the optimization of a cover-like structure with the directions of the pressure loads perpendicular to the loading surface was considered. Figure 11 Figure 11. Design domain and boundary conditions of the cover-like structure. Figure 11. Design domain and boundary conditions of the cover-like structure.
The optimized topology for the cover-like structure is shown in Figure 12, and the final objective function is C = 132.68 Nmm. The elemental relative density of the triangle element in the design domain is shown in Figure 13. As seen in Figure 12, the optimized topology for the cover-like structure is an arched structure, which is similar to the optimized solution obtained by Zhang [12] and Neofytou [23]. Therefore, it can be seen that the proposed method is effective in solving the topology optimization problem with designdependent loads.
Symmetry 2021, 13, 1976 13 of 19 The optimized topology for the cover-like structure is shown in Figure 12, and the final objective function is 132.68 Nmm C = . The elemental relative density of the triangle element in the design domain is shown in Figure 13. As seen in Figure 12, the optimized topology for the cover-like structure is an arched structure, which is similar to the optimized solution obtained by Zhang [12] and Neofytou [23]. Therefore, it can be seen that the proposed method is effective in solving the topology optimization problem with design-dependent loads.

Example 3
The design domain and boundary conditions of the structure in this example are similar to those in Section 5.1, while the position and directions of the initial pressure loads are different. As shown in Figure 14, the bridge-like structure was subjected to constant initial pressure loads =10 N/mm P on the top, left, and right edges and the pressure load direction was perpendicular to the loading surface. The final volume fraction was set to 50%.  The optimized topology for the cover-like structure is shown in Figure 12, and the final objective function is 132.68 Nmm C = . The elemental relative density of the triangle element in the design domain is shown in Figure 13. As seen in Figure 12, the optimized topology for the cover-like structure is an arched structure, which is similar to the optimized solution obtained by Zhang [12] and Neofytou [23]. Therefore, it can be seen that the proposed method is effective in solving the topology optimization problem with design-dependent loads.

Example 3
The design domain and boundary conditions of the structure in this example are similar to those in Section 5.1, while the position and directions of the initial pressure loads are different. As shown in Figure 14, the bridge-like structure was subjected to constant initial pressure loads =10 N/mm P on the top, left, and right edges and the pressure load direction was perpendicular to the loading surface. The final volume fraction was set to 50%.

Example 3
The design domain and boundary conditions of the structure in this example are similar to those in Section 5.1, while the position and directions of the initial pressure loads are different. As shown in Figure 14, the bridge-like structure was subjected to constant initial pressure loads P = 10 N/mm on the top, left, and right edges and the pressure load direction was perpendicular to the loading surface. The final volume fraction was set to 50%.
Symmetry 2021, 13,1976 13 of 19 The optimized topology for the cover-like structure is shown in Figure 12, and the final objective function is 132.68 Nmm C = . The elemental relative density of the triangle element in the design domain is shown in Figure 13. As seen in Figure 12, the optimized topology for the cover-like structure is an arched structure, which is similar to the optimized solution obtained by Zhang [12] and Neofytou [23]. Therefore, it can be seen that the proposed method is effective in solving the topology optimization problem with design-dependent loads.

Example 3
The design domain and boundary conditions of the structure in this example are similar to those in Section 5.1, while the position and directions of the initial pressure loads are different. As shown in Figure 14, the bridge-like structure was subjected to constant initial pressure loads =10 N/mm P on the top, left, and right edges and the pressure load direction was perpendicular to the loading surface. The final volume fraction was set to 50%.   The bridge-like structure optimum solution and the elemental relative density of the design domain are shown in Figures 15 and 16, and the final compliance of this optimum topology is C = 160.88 Nmm. As seen in Figure 15, the curvature center of the curved loading surface is below the structure. As the spherical shape is theoretically an ideal structure for a pressure vessel, this optimized structure is also intuitively the expected solution.
Symmetry 2021, 13,1976 14 of 19 The bridge-like structure optimum solution and the elemental relative density of the design domain are shown in Figures 15 and 16, and the final compliance of this optimum topology is 160.88 Nmm C = . As seen in Figure 15, the curvature center of the curved loading surface is below the structure. As the spherical shape is theoretically an ideal structure for a pressure vessel, this optimized structure is also intuitively the expected solution.

Example 4
This piston head structure was introduced in literature [21], as shown in Figure 17.
The design domain is a rectangle whose size 120 mm 40 mm × , and it was subjected to an initial pressure of =10 N/mm P applied on its top edge. The direction of the pressure loads was perpendicular to the loading surface. In Figure 17, the degrees of freedom along the horizontal direction are constrained on the left and right sides, and the central point of the bottom edge of the structure is fully constrained. A volume fraction of 40% was applied.  The bridge-like structure optimum solution and the elemental relative density of the design domain are shown in Figures 15 and 16, and the final compliance of this optimum topology is 160.88 Nmm C = . As seen in Figure 15, the curvature center of the curved loading surface is below the structure. As the spherical shape is theoretically an ideal structure for a pressure vessel, this optimized structure is also intuitively the expected solution.

Example 4
This piston head structure was introduced in literature [21], as shown in Figure 17.
The design domain is a rectangle whose size 120 mm 40 mm × , and it was subjected to an initial pressure of =10 N/mm P applied on its top edge. The direction of the pressure loads was perpendicular to the loading surface. In Figure 17, the degrees of freedom along the horizontal direction are constrained on the left and right sides, and the central point of the bottom edge of the structure is fully constrained. A volume fraction of 40% was applied.

Example 4
This piston head structure was introduced in literature [21], as shown in Figure 17. The design domain is a rectangle whose size 120 mm × 40 mm, and it was subjected to an initial pressure of P = 10 N/mm applied on its top edge. The direction of the pressure loads was perpendicular to the loading surface. In Figure 17, the degrees of freedom along the horizontal direction are constrained on the left and right sides, and the central point of the bottom edge of the structure is fully constrained. A volume fraction of 40% was applied.
Symmetry 2021, 13,1976 14 of 19 The bridge-like structure optimum solution and the elemental relative density of the design domain are shown in Figures 15 and 16, and the final compliance of this optimum topology is 160.88 Nmm C = . As seen in Figure 15, the curvature center of the curved loading surface is below the structure. As the spherical shape is theoretically an ideal structure for a pressure vessel, this optimized structure is also intuitively the expected solution.

Example 4
This piston head structure was introduced in literature [21], as shown in Figure 17.
The design domain is a rectangle whose size 120 mm 40 mm × , and it was subjected to an initial pressure of =10 N/mm P applied on its top edge. The direction of the pressure loads was perpendicular to the loading surface. In Figure 17, the degrees of freedom along the horizontal direction are constrained on the left and right sides, and the central point of the bottom edge of the structure is fully constrained. A volume fraction of 40% was applied.   Figure 18 shows the structural topology of piston head structure subjected to designdependent loads by using the proposed method, where the compliance value of this optimized design is C = 160.99 Nmm. As shown in Figure 18, from an overall look, the gradient-free optimized solution in this study is similar to the design from Emmendoerfer et al. [21]. In addition, the triangle element grid of the design domain and the elemental relative density are shown in Figure 19.
Symmetry 2021, 13,1976 15 of 19 Figure 18 shows the structural topology of piston head structure subjected to designdependent loads by using the proposed method, where the compliance value of this optimized design is 160.88 Nmm C = . As shown in Figure 18, from an overall look, the gradient-free optimized solution in this study is similar to the design from Emmendoerfer et al. [21]. In addition, the triangle element grid of the design domain and the elemental relative density are shown in Figure 19.

Example 5
The last example considered a two-point supported beam structure, whose sides were 120 mm 30 mm × , where the friction contact effect and geometrical nonlinearity were considered. Generally speaking, for large deformation contact problems, the contact position of a structure changes when a load is applied. Therefore, to ensure the consistency of the structure in the initial position, we added a non-design domain part, whose sides were 120 mm 1 mm × , at the bottom of the structure, as shown in Figure 20. The beam was placed on two semi-circular supports with a radius of An initial pressure P was applied at the top edge of the structure, and the directions of the pressure loads were fixed. The allowable material volume fraction was 50%.  Symmetry 2021, 13,1976 15 of 19 Figure 18 shows the structural topology of piston head structure subjected to designdependent loads by using the proposed method, where the compliance value of this optimized design is 160.88 Nmm C = . As shown in Figure 18, from an overall look, the gradient-free optimized solution in this study is similar to the design from Emmendoerfer et al. [21]. In addition, the triangle element grid of the design domain and the elemental relative density are shown in Figure 19.

Example 5
The last example considered a two-point supported beam structure, whose sides were 120 mm 30 mm × , where the friction contact effect and geometrical nonlinearity were considered. Generally speaking, for large deformation contact problems, the contact position of a structure changes when a load is applied. Therefore, to ensure the consistency of the structure in the initial position, we added a non-design domain part, whose sides were 120 mm 1 mm × , at the bottom of the structure, as shown in Figure 20. The beam was placed on two semi-circular supports with a radius of An initial pressure P was applied at the top edge of the structure, and the directions of the pressure loads were fixed. The allowable material volume fraction was 50%.

Example 5
The last example considered a two-point supported beam structure, whose sides were 120 mm × 30 mm, where the friction contact effect and geometrical nonlinearity were considered. Generally speaking, for large deformation contact problems, the contact position of a structure changes when a load is applied. Therefore, to ensure the consistency of the structure in the initial position, we added a non-design domain part, whose sides were 120 mm × 1 mm, at the bottom of the structure, as shown in Figure 20. The beam was placed on two semi-circular supports with a radius of r = 30 mm. Both supports were rigid and had a friction coefficient of f c = 0.15 with the surface of the beam. An initial pressure P was applied at the top edge of the structure, and the directions of the pressure loads were fixed. The allowable material volume fraction was 50%.
Symmetry 2021, 13,1976 15 of 19 Figure 18 shows the structural topology of piston head structure subjected to designdependent loads by using the proposed method, where the compliance value of this optimized design is 160.88 Nmm C = . As shown in Figure 18, from an overall look, the gradient-free optimized solution in this study is similar to the design from Emmendoerfer et al. [21]. In addition, the triangle element grid of the design domain and the elemental relative density are shown in Figure 19.

Example 5
The last example considered a two-point supported beam structure, whose sides were 120 mm 30 mm × , where the friction contact effect and geometrical nonlinearity were considered. Generally speaking, for large deformation contact problems, the contact position of a structure changes when a load is applied. Therefore, to ensure the consistency of the structure in the initial position, we added a non-design domain part, whose sides were 120 mm 1 mm × , at the bottom of the structure, as shown in Figure 20. The beam was placed on two semi-circular supports with a radius of An initial pressure P was applied at the top edge of the structure, and the directions of the pressure loads were fixed. The allowable material volume fraction was 50%.  The optimized solutions for the structure subjected to different values of pressure loads, P = 10 N/mm and P= 3.2 KN/mm, are presented in Figure 21a,b, respectively. As seen in Figure 21, the two optimized topologies for the cases with the two loads have distinct differences. The length of the solid material at the bottom of the beam structure in Figure 21a is significantly greater than that of the structure in Figure 21b. The reason for this is that with the increase in the pressure load, the contact point between the beam and the two semi-circular supports changes inwards.  Figure 21a,b, respectively. As seen in Figure 21, the two optimized topologies for the cases with the two loads have distinct differences. The length of the solid material at the bottom of the beam structure in Figure 21a is significantly greater than that of the structure in Figure 21b. The reason for this is that with the increase in the pressure load, the contact point between the beam and the two semi-circular supports changes inwards.  The deformed configurations of the optimal solutions in Figure 21 for the two cases are shown in Figure 22a,b. The deformation of the structure for =3.2 KN/mm P was significantly greater than that of the structure for =10 N/mm P ( max =3.38 mm u vs. max =0.027 mm u ). This means that for the topology optimization with design-dependent loads, the large pressure loads have an important influence on the optimal material layout. Here, it should be noted that for the geometrical nonlinearity and contact problems considering the friction coefficient, it would be difficult to analytically derive the sensitivity of the objective function. According to the optimized solutions in Figure 21, this reveals that it is also effective to implement the gradient-free method to deal with the complex topology optimization problem with design-dependent loads.
(a) The deformed configurations of the optimal solutions in Figure 21 for the two cases are shown in Figure 22a,b. The deformation of the structure for P= 3.2 KN/mm was significantly greater than that of the structure for P = 10 N/mm (u max = 3.38 mm vs. u max = 0.027 mm). This means that for the topology optimization with design-dependent loads, the large pressure loads have an important influence on the optimal material layout. Here, it should be noted that for the geometrical nonlinearity and contact problems considering the friction coefficient, it would be difficult to analytically derive the sensitivity of the objective function. According to the optimized solutions in Figure 21, this reveals that it is also effective to implement the gradient-free method to deal with the complex topology optimization problem with design-dependent loads. max loads, the large pressure loads have an important influence on the optimal material layout. Here, it should be noted that for the geometrical nonlinearity and contact problems considering the friction coefficient, it would be difficult to analytically derive the sensitivity of the objective function. According to the optimized solutions in Figure 21, this reveals that it is also effective to implement the gradient-free method to deal with the complex topology optimization problem with design-dependent loads.

Conclusions
In this paper, an effective gradient-free topology optimization strategy is adopted to deal with continuum structures subjected to design-dependent loads. Based on the MFSE model, it is convenient to generate a body-fitted finite element mesh and to identify the loading surface for the application of the design-dependent boundary loads. Furthermore, the design variables used to describe the structure topology are significantly reduced. The topology optimization problem is solved by using the sequential Kriging-based optimization algorithm with a self-adaptive design space adjusting strategy. Numerical examples show the validity and effectiveness of the proposed method.
Although more than 2000 finite element analyses are required to complete the proposed gradient-free topology optimization, it will have developing prospect in practical engineering since the proposed method does not need sensitivity information. With the rapid development of computer hardware and parallel processing technology, the computational costs of the proposed method will not be a heavy burden shortly. It can be easily implemented to solve complex nonlinear problems considering design-dependent loads, contact effects, and geometrical nonlinearity.

Conflicts of Interest:
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Conclusions
In this paper, an effective gradient-free topology optimization strategy is adopted to deal with continuum structures subjected to design-dependent loads. Based on the MFSE model, it is convenient to generate a body-fitted finite element mesh and to identify the loading surface for the application of the design-dependent boundary loads. Furthermore, the design variables used to describe the structure topology are significantly reduced. The topology optimization problem is solved by using the sequential Kriging-based optimization algorithm with a self-adaptive design space adjusting strategy. Numerical examples show the validity and effectiveness of the proposed method.
Although more than 2000 finite element analyses are required to complete the proposed gradient-free topology optimization, it will have developing prospect in practical engineering since the proposed method does not need sensitivity information. With the rapid development of computer hardware and parallel processing technology, the computational costs of the proposed method will not be a heavy burden shortly. It can be easily implemented to solve complex nonlinear problems considering design-dependent loads, contact effects, and geometrical nonlinearity.