1. Introduction
Multiscale analysis has become an important method for the study of materials with microscale structures, such as polymer composites and porous materials with microscale structures [
1]. FE
2 is a two-scale method used to simulate and describe the macroscale structural behavior of material microstructure [
2]. It is a two-layer nested finite element analysis method that uses finite element calculation at macroscale and microscale. The difference is that FE
2 only establishes constitutive equations on the microstructure [
3]. In the multiscale FE
2 model, the periodic boundary conditions are imposed on the representative volume element (RVE) at the Gaussian integration point of each macroscale element. The Newton–Raphson method is used to update the displacement field of the overall structure [
4,
5,
6,
7]. In the nearly three decades since its proposal, the FE
2 method has undergone numerous developments and applications, from the dimensional design of the RVE to interlayer scale transitions [
8,
9] and from quasi-static loading and mechanical coupling problems to multi-physics field phenomena [
10,
11]. Feyel [
12] applied FE
2 to the study of the viscoelastic behavior of long fibers SiC/Ti. Tikarrouchine [
13] solved the problem of the thermal–mechanical coupling of short glass fibers based on FE
2. Kaczmarczyk and Bacigalupo [
14,
15,
16] developed higher-order homogenization schemes.
Traditional FE
2 is very time-consuming because it requires iterative solution. At the same time, programming of the model’s constitutive law is also required. Tan [
17] proposed an improved Direct FE
2 model based on the FE
2 method, which uses multi-point constraints (MPCs) to establish equations between microscale RVEs and macroscale elements. It makes the internal virtual work for microscale FE analysis equal to that of macroscale FE analysis by scaling the sizes of RVEs at Gaussian integration points [
18]. Direct FE
2 eliminates the two-layer nesting problem of the traditional FE
2 method and greatly reduces the data transfer cost, computational cost, and memory cost. In addition, the method can be directly performed on existing mature commercial FE codes (e.g., ABAQUS) without extensive programming knowledge and constitutive modeling expertise, which facilitates the use of FE
2 and expands its application [
19]. The method also has good applications in geometric nonlinear and material nonlinear models [
20,
21,
22].
Topology-optimized design has made significant progress in the last two decades, and the demand for materials with space-varying properties is increasing. Traditional topology optimization methods, such as the variable density method [
23], the evolutionary structural optimization (ESO) method [
24], and the level set method [
25], focus on a single scale. The variable density method, a popular and classical topology optimization method based on the homogenization method, varies element stiffness using a density function. It includes two methods: the solid isotropic material with penalization model (SIMP) [
26] and the rational approximation of material properties model (RAMP) [
27]. The ESO method is a sensitivity topology optimization method based on the von Mises stress criterion that gradually removes inefficient and ineffective material regions from the structure. Xie [
28] proposed the bidirectional evolutionary structural optimization (BESO) method, which adds materials based on previous optimization results. The level set method uses a level set function to determine the initial material distribution and calculates the cell sensitivity using the Lagrangian function [
29].
Multiscale structural topology optimization is a promising technique that enables the design of materials at the microscale instead of just the macroscale. In recent years, concurrent multiscale topological optimization design has become increasingly popular among research scholars. One example of this approach is the concurrent multiscale optimization design proposed by Xia [
30], which optimizes the structure based on the FE
2 method and using the BESO optimization method. Rodrigues [
31] proposed a hierarchical algorithm to optimize material distribution, including local microstructure. Coelho [
32] extended this algorithm to 3D structure optimization. In order to expand the degree of design freedom and propose innovative structures, some scholars have designed lightweight single-cell structures with gradients [
33], heterogeneous periodic structures [
34] and three-scale structures [
35]. Although many scholars have proposed different topology optimization methods, the low computational efficiency and the difficulty of numerical implementation are still problems that need to be solved urgently.
To further improve the efficiency of multiscale optimization, different techniques can be employed. One such technique is the Direct FE2 algorithm, which simplifies the two FE calculations at microscale and macroscale into a single microscale FE calculation. This can greatly improve computational efficiency, but there is a potential issue of checkerboard patterns arising during concurrent multiscale optimization using this method. The main purpose of this paper is to solve the checkerboard problem in the multiscale topology optimization based on Direct FE2, and reconstruct the results to provide reference for structural design. This is of great significance for the practical engineering application of multiscale topology optimization based on Direct FE2.
In this article, we have employed filtering techniques to address the checkerboard pattern issue that may arise when using the Direct FE
2 method for concurrent multiscale optimization. Filtering is a commonly used method in topology optimization. Jantos [
36] used filtering to smooth fiber pathways when optimizing anisotropic materials. In parallel algorithms for topology optimization, a filter based on partial differential equations is more efficient than a density filter [
37]. Yang [
38] proposed a massively efficient filter for topology optimization utilizing the splitting of tensor product structure. Filtering methods are also commonly employed in topology optimization for metal additive manufacturing [
39]. Additionally, we have utilized an interpolation-based approach to reconstruct the optimized results. The proposed iso-position filtering improves the filtering efficiency, and the reconstructed model can serve as a reference for structural design. The multiscale topology optimization based on Direct FE
2 can be implemented on finite element commercial software, thus significantly reducing the professional expertise required for multiscale topology optimization and facilitating its widespread application.
This paper is organized as follows:
Section 2 introduces concurrent multiscale topology optimization, including the multiscale topology optimization based on Direct FE
2.
Section 3 proposes the filtering and reconstruction method to reconstruct the optimized results. The proposed method is then verified by three numerical examples in
Section 4. Conclusions are made in
Section 5.
3. Reconstruction of Concurrent Multiscale Optimization Based on Direct FE2
The concurrent multiscale topology optimization based on Direct FE
2 may result in a checkerboard pattern. At the same time, the optimization result of Direct FE
2 contains only the densities of RVE located at the Gauss points of macroscale elements, which are not necessarily continuous with each other and represent only the tendency of the topological variations. For manufacturing purposes, denser material points need to be considered. Some image-based interpolation schemes such as [
41] are developed for generating intermediate microstructures. Here, a reconstruction method based on the interpolation of element densities is proposed, which provides a direct way to show the reconstruction result of topology optimization.
The reconstruction scheme is shown in
Figure 1. After optimization of the Direct FE
2 model, the density results
are read from the result file (for Abaqus, the .odb file). Then, all element densities for each macroscale element are interpolated: see
Figure 2. For rectangular elements with
Gauss quadrature points, the following bilinear interpolation function is used.
where
is the interpolated densities and
are the center coordinates of the
RVEs whose densities are interpolated. The coefficients
,
,
, and
can be solved according to the densities of RVE located at the Gauss point of macroscale elements.
In order to avoid the checkerboard problem and mesh dependency [
42], the density filter of elements is performed after interpolating element densities. Different from the traditional density filter method, which takes into account all the elements that fall within the filter radius, an iso-position filter method is proposed here to speed up the filter procedure. There are two steps for the iso-position filter: see
Figure 3. Firstly, the RVEs within the filter radius
are searched according to the center coordinates of the RVEs. Secondly, the filtered densities of each element in the RVE are calculated using the following formula.
where
is the filtered density and
is the interpolated density of the element located in the same position of the RVE with element
.
is the number of iso-position elements whose distance with element
is less than the filter radius
. The weight
is defined by
where
is the distance of element
and
.
Finally, the filtered densities are sorted and the reconstruction model is built based on ordered element densities and volume targets. It should be noted that the interpolated densities may be negative values when the densities of the four RVEs vary widely because extrapolation is used here. However, this does not affect the reconstruction result. It just represents a rather lower density than the positive densities. These elements will be ranked in a later place than those with positive densities after the sorting of the densities.
The pseudo-code of reconstruction is listed in Algorithm 1. There are two main loops in this algorithm. Line 2 to line 6 is the first main loop, which interpolates the densities of all elements. The coefficients
,
,
, and
are the same for the
elements located in the same position of the RVE in one macroscale element. Therefore, the coefficients only need to be solved once for each RVE element, and the loop of line 3 is outside the loop of line 5. The second main loop, which filters the element densities, is from line 7 to line 16. The weight
is the same for the iso-position elements of the RVE
and
. Therefore, the loop of line 10 is outside the loop of line 13. In order to reduce the computational cost, the neighbor RVEs are searched in advance in line 8. Thus, the weight formula Equation (22) is directly replaced by
in line 11.
Algorithm 1 Reconstruction of optimization result of Direct FE2 |
1: | Read density from the optimization result file |
2: | For each element in the macroscale mesh: |
3: | For each element in the RVE mesh: |
4: | Calculate the interpolating coefficients , , , and |
5: | For RVEs of the macroscale element: |
6: | Calculate the interpolated densities |
7: | For each interpolated RVE : |
8: | Find the neighbor RVEs whose distance with RVE is less than filter radius |
9: | Initialize density and the sum of weight |
10: | For each neighbor RVE : |
11: | Calculate the weight |
12: | Update the sum of weight |
13: | For each element in RVE : |
14: | Update density by iso-position element of RVE : |
15: | For each element in RVE : |
16: | Divide the density by the sum of weight |
17: | Sort the filtered densities |
18: | Build the reconstruction model based on ordered element densities and volume target |
4. Numerical Examples
In this section, some 2D numerical examples are presented to demonstrate the capabilities of topology optimization based on Direct FE2. The commercial software, Abaqus 2020, is used with CPS4 plane stress elements for all of these simulations. The built-in SIMP method is used for the topology optimization simulation. The target volume fraction for the topology optimization is and the penalization power is set to . The default values are used for other algorithm parameters. The stopping criterion in these examples is that both the change of density between the two iterations was less than and the change of objective function between the two iterations was less than . Moreover, the maximum number of iterations is .
The size of the macroscale elements is
, and it is assumed that each macroscale element contains
RVEs. The thickness of the multiscale structure is
. Thus, the thickness of the RVE mesh is set to
according to Equation (8). Due to the applied periodic boundary conditions, an initial guess design of the RVE has to be defined to avoid a uniformly distributed sensitivity field [
43]. Here, the initial design is one hole in each RVE, and the radius of the hole is a quarter of the edge length of the RVE. The RVE is discretized into a structured mesh with
elements. The Young’s modulus and Poisson’s ratio of the cellular material are
and
, respectively. Due to the macro elements not contributing to the stiffness in the Direct FE
2 algorithm, Young’s modulus of macroscale element is set to
. The results of Direct FE
2 are compared with the direct numerical simulation (DNS) model.
4.1. Bridge-Type Structure Model
In this first example, the bridge-type structure in
Figure 4 is studied, where a concentrated force is loaded at the middle point of the bottom side. The length of the whole structure
is
, and the height
is
. The Direct FE
2 model is discretized into
macroscale elements, and each contains
RVEs located at the Gauss points. Each RVE has
elements, as shown in
Figure 5.
The contour plot of densities after optimization is shown in
Figure 6. To better show the results, only the elements with the top
of the density are drawn in the whole model picture. Four typical macroscale elements are zoomed in and plotted at the bottom of the whole model picture. The
,
,
, and
in the figure represent the retained volume fraction of each RVE.
It can be seen that the volume fractions of each RVE are different from one other in
Figure 6. Therefore, the Direct FE
2 has more freedom in the local optimization and is more flexible than other concurrent topology optimization methods [
30]. However, due to the minimum density being set to
to prevent the model from a singularity, the RVEs with volume fraction of
should be removed totally according to the optimization result. This will lead to the checkerboard problem or “gap” problem, which will be shown in the following explanation. Therefore, the density filter in
Section 3 is necessary. The density filter will “fill in” the “gap” and the reconstruction result are shown in
Figure 7. Some local details of the reconstructed Bridge-type model are also plotted, and it can be seen that the local details are consistent with the force flow of the structure.
The strain energy in the iterative process of the Direct FE
2 model and the DNS model is shown in
Figure 8. The target volume fraction of the Direct FE
2 model and the DNS model for the topology optimization is
. The strain energy of the Direct FE
2 model and the DNS model before optimization is
and
, respectively. After optimization, the strain energy of the Direct FE
2 model and the DNS model becomes
and
, respectively. Due to the initial strain energies of the Direct FE
2 model and the DNS model being different, the strain energy is normalized by the initial strain energy of each computational model [
17]. It shows that the strain energy of the two models decreases steadily with the increase in the number of iterations. It can be seen from the figure that the Direct FE
2 model converges faster and that the normalized strain energy of the two models is the same (
) after optimization.
The topology optimization results of Direct FE
2 and DNS are shown in
Figure 9. The optimized pattern of DNS contains some fibrous structure. This is due to the fact that the DNS model has too many microscale elements and the optimization results are heavily mesh-dependent.
It should be noted that the Direct FE
2 simulation required only about
min using
processors to complete, whereas the DNS model took about
h
min with
processors to reach optimization completion. The number of elements, the number of iterations, and the calculation times of Direct FE
2 and DNS are shown in
Table 2.
Four representative filter radiuses,
are considered to illustrate the choice of filter radius. The results of the four test cases are shown in
Figure 10a–d, correspondingly. It shows that the density filter is very helpful for the reconstruction of the optimization result of the Direct FE
2 model. The “gap” problem occurs when the filter radius is too small, and the boundaries of the reconstruction result will be rough. The boundaries will become smoother as the filter radius increases. When the filter radius is too large, more macroscale materials are removed because the densities of microscale elements tend to be the same. In this example, a filter radius between
and
is recommended. The topological evolution of the Direct FE
2 model with filter radius
is shown in
Figure 11.
4.2. MBB Beam Model
A classical topology optimization example, the Messerschmitt–Bolkow–Blohm (MBB) beam structure, is considered in this example. The height of the whole structure
is set to
, and the length
is set to
. A concentrated force of
is loaded at the middle point of the top side. Due to the symmetry of the MBB beam, only one half of the beam is analyzed, as shown in
Figure 12. For the Direct FE
2 model, the half of the MBB beam is discretized into
macro elements. Each macroscale element has
RVEs at the Gaussian integration point and each RVE is discretized into a structured mesh with
elements. The volume fraction for the topology optimization of the design domain is defined as
.
The contour plot of densities after optimization is shown in
Figure 13. The elements with the top
of the density are drawn in the whole model picture. Four typical macroscale elements are zoomed in and plotted at the bottom of the whole model picture. The
,
,
, and
in the figure represent the retained volume fraction of each RVE.
The reconstruction result with filter radius
is shown in
Figure 14. The optimization result shows that uniaxial materials may be sufficient at the major branches of the structure, while anisotropic materials have to be used at the joints of the major branches in order to provide the higher structural stiffness required due to the more complex load conditions.
The strain energy in the iterative process of the Direct FE
2 model and the DNS model is shown in
Figure 15. The target volume fraction of the Direct FE
2 model and the DNS model for the topology optimization is
. The strain energy of the Direct FE
2 model and the DNS model before optimization is
and
, respectively. After optimization, the strain energy of the Direct FE
2 model and the DNS model becomes
and
, respectively. Due to the initial strain energies of the Direct FE
2 model and the DNS model being different, the strain energy is normalized by the initial strain energy of each computational model. The strain energy decreases steadily with the increase in the number of iterations. The Direct FE
2 model converges faster than the DNS model, and the Direct FE
2 model has relatively lower strain energy compared with the DNS model (normalized strain energy of
and
, respectively).
The topology optimization results of Direct FE
2 and DNS are shown in
Figure 16. The optimized pattern of DNS contains a more fibrous structure because the DNS model has too many microscale elements and the optimization results are heavily mesh-dependent.
Direct FE
2 simulation required only about
min using
processors to complete the optimization, whereas the DNS model took about
h
min with
processors to reach optimization completion. The number of elements, the number of iterations, and the calculation times of Direct FE
2 and DNS are shown in
Table 3.
Four representative filter radiuses,
are considered to illustrate the choice of filter radius: see
Figure 17. A filter radius between
and
is recommended for the MBB beam model. The topological evolution of the Direct FE
2 model with filter radius
is shown in
Figure 18.
4.3. Loaded Knee Model
To further demonstrate the optimization capability of the Direct FE
2 method, we considered the loaded knee model. The dimensional details and schematic configuration of the
-shaped structure are indicated in
Figure 19. The upper end of the loaded knee structure is fully restrained, and the concentrated force is loaded at the midpoint of the right foot side of the loaded knee structure. The geometry of dimensions
is set to
, and the concentrated force
is set to
. The volume fraction of the design domain is defined as
.
The contour plot of densities after optimization is shown in
Figure 20. To better show the results, only the elements with the top
of the density are drawn in the whole model picture. Four typical macroscale elements are zoomed in and plotted at the bottom of the whole model picture. The
,
,
, and
in the figure represent the retained volume fraction of each RVE.
Figure 20 further verifies the flexibility of the optimization method. The volume fractions of the RVEs are different from one other. The reconstruction model with some local details is plotted in
Figure 21, and it can be seen that the local details are consistent with the force flow of the structure.
The target volume fraction of the Direct FE
2 model and the DNS model for the topology optimization is
. The strain energy of the Direct FE
2 model and DNS model before optimization is
and
, respectively. After optimization, the strain energy of the Direct FE
2 model and DNS model becomes
and
, respectively. Due to the initial strain energies of the Direct FE
2 model and the DNS model being different, the strain energy is normalized by the initial strain energy of each computational model. The normalized strain energy in the iterative process of the Direct FE
2 model and the DNS model is shown in
Figure 22. The normalized strain energy of the two models is
and
after optimization.
The topology optimization results of Direct FE
2 and DNS are shown in
Figure 23. The optimized pattern of DNS contains a more fibrous structure because the DNS model has too many microscale elements and the optimization results are heavily mesh-dependent.
Direct FE
2 simulation required only about
min using
processors to complete the optimization, whereas the DNS model took about
h with
processors to reach optimization completion. The number of elements, number of iterations and calculation times of Direct FE
2 and DNS are shown in
Table 4.
Four representative filter radiuses,
are considered to illustrate the choice of filter radius: see
Figure 24. The “gap” problem occurs when the filter radius is too small. A filter radius between
and
is recommended for the loaded knee model. The topological evolution of the Direct FE
2 model with filter radius
is shown in
Figure 25.
5. Conclusions
The main purpose of this paper is to address the issues that arise in practical applications of multiscale topology optimization. Traditional multiscale topology optimization methods are generally time-consuming and difficult to implement in programming. However, the concurrent multiscale topology optimization based on Direct FE2 can greatly reduce the computational cost and can be implemented on general commercial FE software, greatly reducing the threshold for multiscale topology optimization. Nevertheless, the concurrent multiscale topology optimization based on Direct FE2 can encounter the checkerboard problem. In this paper, the checkerboard problem is addressed, and the model is reconstructed to make the Direct FE2-based multiscale topology optimization method more practical.
In order to solve the checkerboard problem in the results of Direct FE2-based multiscale topology optimization, we reconstructed the results using bilinear interpolation and filtering. By extrapolating the density values of the RVEs at four integration points, the density values of all RVEs within the macroscale element can be obtained. Then, the interpolated density is filtered by performing a weighted average of the element density within a circular area with a radius of . This is the key step to eliminate the checkerboard problem. In order to improve the efficiency of filtering, this paper proposes an iso-position filtering method. When filtering the density of a certain element, only the weighted average of the element density at the same position in the RVE is performed, instead of all the elements within the circular area. Finally, the reconstructed structure is determined according to the filtered density of the elements.
Some numerical examples are also presented to verify the reconstruction method proposed in this paper. The multiscale topology optimization based on Direct FE2 converges faster and yields consistent normalized strain energy after optimization, in comparison with DNS. From the numerical examples, it can be seen that a filter radius between 40 cm and 60 cm is appropriate. The size of the macroscale elements is 20 cm. Therefore, the recommended filter radius is two to three times the size of the macroscale elements. It should be noted that in addition to the size of the macroscale elements, the overall size of the structure also needs to be considered when selecting the filtering radius. When the filtering radius is too small, the checkerboard pattern problem still exists. When the filtering radius is too large, the density of each element tends to be uniform, making it difficult to effectively identify the elements that need to be removed.
Overall, this study proposes a novel algorithm for concurrent multiscale topology optimization based on Direct FE2 that addresses the checkerboard problem and significantly reduces the computational cost. The proposed algorithm can be easily implemented on general commercial FE software, making it accessible to a broader range of researchers and engineers. Future research could focus on applying the proposed algorithm to more complex problems and exploring its potential for practical engineering applications.