Next Article in Journal
Oscillation Criteria for Qusilinear Even-Order Differential Equations
Next Article in Special Issue
Two-Dimensional Equivalent Models in the Analysis of a Multibody Elastic System Using the Finite Element Analysis
Previous Article in Journal
Self-Similarity Principle and the General Theory of Fractal Elements: How to Fit a Random Curve with a Clearly Expressed Trend?
Previous Article in Special Issue
Stability Analysis of the Rational Solutions, Periodic Cross-Rational Solutions, Rational Kink Cross-Solutions, and Homoclinic Breather Solutions to the KdV Dynamical Equation with Constant Coefficients and Their Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Reconstruction Approach for Concurrent Multiscale Topology Optimization Based on Direct FE2 Method

1
School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, China
2
Department of Mechanical Engineering, National University of Singapore, Singapore 119260, Singapore
3
International Machinery Center, Department of Mechanical Engineering, Xi’an Jiaotong University, Xi’an 710049, China
4
Shanghai Aerospace Control Technology Institute, Shanghai 201109, China
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(12), 2779; https://doi.org/10.3390/math11122779
Submission received: 31 May 2023 / Revised: 17 June 2023 / Accepted: 18 June 2023 / Published: 20 June 2023
(This article belongs to the Special Issue Applied Mathematics and Continuum Mechanics)

Abstract

:
The rapid development of material science is increasing the demand for the multiscale design of materials. The concurrent multiscale topology optimization based on the Direct FE2 method can greatly improve computational efficiency, but it may lead to the checkerboard problem. In order to solve the checkerboard problem and reconstruct the results of the Direct FE2 model, this paper proposes a filtering-based reconstruction method. This solution is of great significance for the practical application of multiscale topology optimization, as it not only solves the checkerboard problem but also provides the optimized full model based on interpolation. The filtering method effectively eliminates the checkerboard pattern in the results by smoothing the element densities. The reconstruction method restores the smoothness of the optimized structure by interpolating between the filtered densities. This method is highly effective in solving the checkerboard problem, as demonstrated in our numerical examples. The results show that the proposed algorithm produces feasible and stable results.

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]. FE2 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 FE2 only establishes constitutive equations on the microstructure [3]. In the multiscale FE2 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 FE2 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 FE2 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 FE2. Kaczmarczyk and Bacigalupo [14,15,16] developed higher-order homogenization schemes.
Traditional FE2 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 FE2 model based on the FE2 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 FE2 eliminates the two-layer nesting problem of the traditional FE2 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 FE2 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 FE2 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 FE2 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 FE2 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 FE2. 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.

2. Concurrent Multiscale Topology Optimization

In this chapter, concurrent multiscale topology optimization based on the FE2 method are reviewed. Moreover, the theory and implementation details of the concurrent multiscale topology optimization based on the Direct FE2 method are discussed.

2.1. The Mathematical Description of Multiscale Topology Optimization

The mathematical description of multiscale topology optimization provided by Theocaris [40] will be adopted in this study. The design variables of the macrostructure and microstructure are denoted by ρ and η , respectively. The mathematical formulation of the multiscale topology optimization problem is as follows.
m a x ρ , η S   m i n u U   1 2 Ω   C i j k l ρ , η ϵ i j u ϵ k l u d Ω l u
where C i j k l is the elastic stiffness tensor and ϵ i j is the strain tensor. l ( u ) is the external force potential energy. S is the design variable space, which can be stated as
S = ρ , η ρ , η = 0   or   1 , Ω   ρ d Ω = V f 1 , Ω   η d Ω = V f 2
where V f 1 and V f 2 are given volume fractions of the macrostructure and microstructure, respectively. The problem 1 can be decomposed into
m a x ρ S 1   m a x η S 2 ρ   m i n u U   1 2 Ω   C i j k l ρ , η ϵ i j u ϵ k l u d Ω l u
S 1 and S 2 are the design variable space for ρ and η respectively. The order of the operators in Equation (3) can be changed to obtain
max ρ S 1   m i n u U   Ω   m a x η S 2   1 2 C i j k l ρ , η ϵ i j u ϵ k l u d Ω l u
The inner maximization problem is designed to find the microstructures with maximum stiffness under macroscale strain, while the outer maximization problem is used to find the macroscale structure with maximum stiffness. The intermediate minimization problem can be solved by a general finite element code.

2.2. Multiscale Topology Optimization Based on FE2

The FE2 method employs nested meshes for the iterative solving of macroscale and microscale structures, where the microscale structures are represented by representative volume elements (RVEs). The constitutive behavior of the macroscale structure is indirectly reflected through the constitutive behavior of the microscale structure, thus eliminating the need for a macroscale constitutive model.
Topology optimization based on the FE2 method extends the original FE2 method by incorporating both macroscale and microscale scale topology optimization. The discrete form of the topology multiscale optimization problem can be expressed as
min U ( ρ , η ) = f T d s . t . K ( ρ , η ) d = f i = 1 N   ρ i V i = V f 1 j = 1 n   η j x v j = V f 2 0 < ρ m i n ρ i 1 , i = 1 , , N 0 < η m i n η j 1 , j = 1 , , n
where U ( ρ , η ) is the structural compliance. K , d , and f are the stiffness matrix, nodal displacements, and external force, respectively. In the FE2 framework, the equation K ( ρ , η ) d = f represents the equilibrium of the macroscale structure, and an iterative method is required to solve it in practical computations. V i and v j , respectively, represent the volume of the macroscale and microscale elements. N and n are the number of macroscale and microscale elements, respectively. The purpose of the minimum value ρ m i n , η m i n is to prevent singularity in the equilibrium problem.
The multiscale topology optimization problem can be decomposed into macroscale and microscale topology optimization according to Equation 4 . The optimization problems for the inner layer Equation 4 constitute a representation of topological optimization at the microscale, while the optimization problems for the outer layer represent topological optimization at the macroscale. Topological optimization at the macroscale can be expressed as
min U ( ρ ) = f T d s . t . K ( ρ , η ) d = f i = 1 N   ρ i V i = V f 1 0 < ρ m i n ρ i 1 , i = 1 , , N
The above equation K ( ρ , η ) d = f implicitly includes the variable η that determines the constitutive relationship of the structure. This equation will be solved through FE2. The optimization of the microstructure can be expressed as follows:
m i n   U i ( η ) = u i T K i η u i s . t . K i η d i = f i   i = 1 n   η i v i = V f 2   0 < η m i n η i 1 , i = 1 , , n
The force f i is calculated based on the macroscale structural strain. The macroscale structure’s strain is applied to the boundary of the RVE through periodic boundary conditions (PBC).
Equations 6 and 7 provide a mathematical description of the discrete form of the topology multiscale optimization problem. It is optimized by embedding it into the FE2 model. As the FE2 method also requires iterative solving, there are three nested loops for multiscale topology optimization based on FE2. This greatly reduces the efficiency of multiscale topology optimization, especially when optimizing large models. Direct FE2 eliminates the iteration between two scales through stiffness scaling, and the macroscale elements are given negligible stiffness, so only the microstructure needs to be optimized. This greatly reduces the computational cost. The following subsection will introduce the multiscale topology optimization based on Direct FE2.

2.3. Multiscale Topology Optimization Based on Direct FE2

The Direct FE2 technique is a highly effective computational approach that integrates two finite element analyses of varying scales into a single one. It employs Representative Volume Elements (RVEs) to depict microstructures and connects the two scale computations using a scale factor of stiffness matrix. By eliminating the iteration between macroscale and microscale simulations, Direct FE2 enhances computational efficiency. Due to the fact that only MPC is needed to establish the connection between macro and micro, it is possible to solve the Direct FE2 model on commonly used commercial finite element software. For a detailed derivation process of Direct FE2, please refer to reference [17]. Here, we will only provide a basic introduction to Direct FE2.
Direct FE2 requires the establishment of two sets of finite element meshes, namely, macroscale and microscale mesh. The macroscale mesh is used to represent the macroscale structure, and a microscale mesh of the RVE is established at the integration points of the macroscale element to simulate the local microscale structure. Direct FE2 assigns a negligible stiffness to the macroscale element and scales the stiffness of the microscale element k ~ i j to conserve the energy of the structure. The scaling factor s of the microscale element stiffness matrix is
s = w α J α 1 V α
where w α is the weight of integration corresponding to Gauss point α . V α represents the volume of the RVE at Gauss point α . J α is the Jacobian determinant for the transformation of physical coordinates to natural coordinates. For the bilinear rectangle elements, J α = V e / 4 , w α = 1 . Therefore, the scale factor for the bilinear rectangle elements is
s = V e 4 V α
where V e is the volume of the macroscale element. After discretization of the virtual work principle, the Direct FE2 equation becomes
K ~ i j δ d ~ i d ~ j = f j δ d j
The symbols with a tilde ‘∼’ indicate the variables related to the RVE. K ~ i j is the scaled stiffness matrix. d ~ j and d j represent the nodal displacements of microscale mesh and macroscale mesh, respectively. f j is the external force acting on macroscale elements. In order to establish the connection between the displacement of macroscale elements and microscale elements, PBC needs to be applied. The displacement of nodes on the RVE boundary is equal to the interpolated displacement within the macroscale elements plus a periodic fluctuation field, which can be stated as
u ~ I = N i ξ , η d i + ϵ ~
Due to the periodicity of the RVE, we then have
ϵ ~ T = ϵ ~ B
where ϵ ~ T and ϵ ~ B are the perturbation of the displacements of nodes at the top edge and bottom edge of RVE, respectively. Then, subtracting the displacement of the top edge from that of the bottom edge will lead to
u ~ T u ~ B = N i ξ T , η T d i N i ξ B , η B d i
where u ~ T and u ~ B represent the displacement of the top and bottom boundaries of the RVE. N i is the shape function of elements. ξ T , η T and ξ B , η B are the natural coordinates of the macroscopic elements corresponding to the upper and lower edges of the RVE, respectively. d i is the nodal displacement of macroscale elements. Similarly, the PBC on the left and right sides of the RVE are expressed as follows.
u ~ R u ~ L = N i ξ R , η R d i N i ξ L , η L d i
After discretization, these PBCs can be expressed as the relationship between macro- and micro-node displacements using the following equation.
d j = L i j d ~ i
By substituting Equation 15 into Equation 10 and eliminating the variation, the equilibrium equation for Direct FE2 can be obtained.
K ~ i j d ~ j = L i j f j
Topology optimization problems can be solved using various methods, and one classic approach is the density-based method. This method represents stiffness as a function of “density”, which enables the optimization of material distribution to achieve optimal performance. One commonly used implementation of the density-based method is the SIMP model, which uses a power law function to relate density to stiffness.
k ~ i j e ρ e = ρ e p k ~ i j e
where k ~ i j e is the stiffness matrix of microscale elements. The index p has been shown to satisfy the following equation in order to meet the Hashin–Shtrikman bounds [23].
p max 2 1 ν , 4 1 + ν   ( for   2D   problems ) , max 15 1 ν 7 5 ν , 3 2 1 ν 1 2 ν   ( for   3D   problems )
where ν is the Poisson’s ratio of the material. Therefore, topology optimization based on Direct FE2 can be expressed as follows.
min   a ( ρ ~ ) = K ~ i j ( ρ ~ ) d ~ i d ~ j = e ρ ~ e p k ~ i j e d ~ i e d ~ j e s . t . K ~ i j ( ρ ~ ) d ~ j = L i j f j V ~ ρ / V ~ 0 = v f 0 < ρ ~ m i n ρ ~ e 1
The topology optimization algorithm in finite element commercial software can be used to solve the above equation. The design domain for multiscale topology optimization based on Direct FE2 is the RVE at the integration point. Table 1 compares different topology optimization models. The DNS model can only perform topology optimization at a single scale, and due to the need to establish a complete computational model, the number of elements is large. The FE2 model requires fewer elements, but due to the need for programming and nested loops, the computational efficiency is low. The Direct FE2 model can use commercial finite element software to perform multiscale topology optimization.

3. Reconstruction of Concurrent Multiscale Optimization Based on Direct FE2

The concurrent multiscale topology optimization based on Direct FE2 may result in a checkerboard pattern. At the same time, the optimization result of Direct FE2 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 FE2 model, the density results ρ 0 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 2 × 2 Gauss quadrature points, the following bilinear interpolation function is used.
ρ ^ = a + b x + c y + d x y
where ρ ^ is the interpolated densities and x , y are the center coordinates of the n × n RVEs whose densities are interpolated. The coefficients a , b , c , and d 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 r f i l t e r 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.
ρ ¯ i = j = 1 N i   w i j ρ ^ j j = 1 N i   w i j
where ρ ¯ i is the filtered density and ρ ^ j is the interpolated density of the element located in the same position of the RVE with element i . N i is the number of iso-position elements whose distance with element i is less than the filter radius r f i l t e r . The weight w i j is defined by
w i j = m a x 0 , r f i l t e r d i , j
where d i , j is the distance of element i and j .
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 a , b , c , and d are the same for the n × n 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 w i j is the same for the iso-position elements of the RVE i and j . 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 w i j = r f i l t e r d i , j in line 11.
Algorithm 1 Reconstruction of optimization result of Direct FE2
1: Read density ρ 0 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 a , b , c , and d
5:       For n × n RVEs of the macroscale element: 
6:         Calculate the interpolated densities ρ ^ = a + b x + c y + d x y
7: For each interpolated RVE i : 
8:   Find the neighbor RVEs whose distance d i , j with RVE i is less than filter radius r f i l t e r
9:   Initialize density ρ ¯ i = 0 and the sum of weight w s u m = 0
10:   For each neighbor RVE j : 
11:     Calculate the weight w i j = r f i l t e r d i , j
12:     Update the sum of weight w s u m = w s u m + w i j
13:     For each element in RVE i : 
14:       Update density ρ ¯ i by iso-position element of RVE j : ρ ¯ i = ρ ¯ i + w i j ρ ^ j
15:   For each element in RVE i : 
16:     Divide the density by the sum of weight ρ ¯ i = ρ ¯ i / w s u m
17: Sort the filtered densities ρ ¯ i
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 50 % and the penalization power is set to 3 . 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 0.005 and the change of objective function between the two iterations was less than 0.001 . Moreover, the maximum number of iterations is 50 .
The size of the macroscale elements is 20   m m × 20   m m , and it is assumed that each macroscale element contains 10 × 10 RVEs. The thickness of the multiscale structure is 5   m m . Thus, the thickness of the RVE mesh is set to 125   m m 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 96 elements. The Young’s modulus and Poisson’s ratio of the cellular material are 70   G P a and 0.3 , respectively. Due to the macro elements not contributing to the stiffness in the Direct FE2 algorithm, Young’s modulus of macroscale element is set to 1 × 10 9   G P a . The results of Direct FE2 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 L is 800   c m , and the height H is 400   c m . The Direct FE2 model is discretized into 40 × 20 macroscale elements, and each contains 2 × 2 RVEs located at the Gauss points. Each RVE has 96 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 50 % 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 V f 1 , V f 2 , V f 3 , and V f 4 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 FE2 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 0.001 to prevent the model from a singularity, the RVEs with volume fraction of 0.10 % 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 FE2 model and the DNS model is shown in Figure 8. The target volume fraction of the Direct FE2 model and the DNS model for the topology optimization is 50 % . The strain energy of the Direct FE2 model and the DNS model before optimization is 1.703 × 10 4   m J and 2.705 × 10 4   m J , respectively. After optimization, the strain energy of the Direct FE2 model and the DNS model becomes 2.795 × 10 5   m J and 4.439 × 10 5   m J , respectively. Due to the initial strain energies of the Direct FE2 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 FE2 model converges faster and that the normalized strain energy of the two models is the same ( = 0.164 ) after optimization.
The topology optimization results of Direct FE2 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 FE2 simulation required only about 49 min using 20 processors to complete, whereas the DNS model took about 19 h 13 min with 20 processors to reach optimization completion. The number of elements, the number of iterations, and the calculation times of Direct FE2 and DNS are shown in Table 2.
Four representative filter radiuses, r f i l t e r = 20   c m ,   40   c m ,   60   c m ,   a n d   80   c m , 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 FE2 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 40   c m and 60   c m is recommended. The topological evolution of the Direct FE2 model with filter radius 60   c m 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 H is set to 400   c m , and the length L is set to 1600   c m . A concentrated force of 1   N 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 FE2 model, the half of the MBB beam is discretized into 40 × 20 macro elements. Each macroscale element has 2 × 2 RVEs at the Gaussian integration point and each RVE is discretized into a structured mesh with 96 elements. The volume fraction for the topology optimization of the design domain is defined as 0.5 .
The contour plot of densities after optimization is shown in Figure 13. The elements with the top 50 % 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 V f 1 , V f 2 , V f 3 , and V f 4 in the figure represent the retained volume fraction of each RVE.
The reconstruction result with filter radius r f i l t e r = 60   c m 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 FE2 model and the DNS model is shown in Figure 15. The target volume fraction of the Direct FE2 model and the DNS model for the topology optimization is 50 % . The strain energy of the Direct FE2 model and the DNS model before optimization is 9.148 × 10 4   m J and 2.783 × 10 4   m J , respectively. After optimization, the strain energy of the Direct FE2 model and the DNS model becomes 1.559 × 10 4   m J and 4.967 × 10 5   m J , respectively. Due to the initial strain energies of the Direct FE2 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 FE2 model converges faster than the DNS model, and the Direct FE2 model has relatively lower strain energy compared with the DNS model (normalized strain energy of 0.170 and 0.178 , respectively).
The topology optimization results of Direct FE2 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 FE2 simulation required only about 41 min using 20 processors to complete the optimization, whereas the DNS model took about 16 h 30 min with 20 processors to reach optimization completion. The number of elements, the number of iterations, and the calculation times of Direct FE2 and DNS are shown in Table 3.
Four representative filter radiuses, r f i l t e r = 20   c m , 40   c m , 60   c m ,   a n d   80   c m , are considered to illustrate the choice of filter radius: see Figure 17. A filter radius between 40   c m and 60   c m is recommended for the MBB beam model. The topological evolution of the Direct FE2 model with filter radius 60   c m is shown in Figure 18.

4.3. Loaded Knee Model

To further demonstrate the optimization capability of the Direct FE2 method, we considered the loaded knee model. The dimensional details and schematic configuration of the L -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 L is set to 800   c m , and the concentrated force F is set to 1   N . The volume fraction of the design domain is defined as 0.5 .
The contour plot of densities after optimization is shown in Figure 20. To better show the results, only the elements with the top 50 % 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 V f 1 , V f 2 , V f 3 , and V f 4 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 FE2 model and the DNS model for the topology optimization is 50 % . The strain energy of the Direct FE2 model and DNS model before optimization is 9.522 × 10 4   m J and 9.903 × 10 4   m J , respectively. After optimization, the strain energy of the Direct FE2 model and DNS model becomes 1.599 × 10 4   m J and 1.747 × 10 4   m J , respectively. Due to the initial strain energies of the Direct FE2 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 FE2 model and the DNS model is shown in Figure 22. The normalized strain energy of the two models is 0.168 and 0.176 after optimization.
The topology optimization results of Direct FE2 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 FE2 simulation required only about 50 min using 20 processors to complete the optimization, whereas the DNS model took about 64 h with 20 processors to reach optimization completion. The number of elements, number of iterations and calculation times of Direct FE2 and DNS are shown in Table 4.
Four representative filter radiuses, r f i l t e r = 20   c m , 40   c m ,   60   c m ,   a n d   80   c m , 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 40   c m and 60   c m is recommended for the loaded knee model. The topological evolution of the Direct FE2 model with filter radius 60   c m 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 r . 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.

Author Contributions

Conceptualization, A.Z., V.B.C.T. and P.L.; methodology, A.Z.; software, A.Z.; validation, A.Z., V.B.C.T. and P.L.; investigation, K.L.; data curation, K.L.; writing—original draft preparation, A.Z.; writing—review and editing, Z.H.; visualization, A.Z. and K.L.; supervision, Z.H.; project administration, Z.H.; funding acquisition, Z.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the China Scholarship Council (CSC), grant number 202006260126.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lakes, R. Materials with structural hierarchy. Nature 1993, 361, 511–515. [Google Scholar] [CrossRef]
  2. Feyel, F. Multiscale FE2 elastoviscoplastic analysis of composite structures. Comput. Mater. Sci. 1999, 16, 344–354. [Google Scholar] [CrossRef]
  3. Feyel, F. A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua. Comput. Methods Appl. Mech. Eng. 2003, 192, 3233–3244. [Google Scholar] [CrossRef]
  4. Ypma, T.J. Historical Development of the Newton-Raphson Method. SIAM Rev. 1995, 37, 531–551. [Google Scholar] [CrossRef] [Green Version]
  5. Smit, R.J.M.; Brekelmans, W.A.M.; Meijer, H.E.H. Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling. Comput. Methods Appl. Mech. Eng. 1998, 155, 181–192. [Google Scholar] [CrossRef]
  6. Hill, R. A Self-Consistent Mechanics of Composite Materials. J. Mech. Phys. Solids 1965, 13, 213–222. [Google Scholar] [CrossRef]
  7. Li, S.; Mao, Y.; Liu, W.; Hou, S. A highly efficient multi-scale approach of locally refined nonlinear analysis for large composite structures. Compos. Struct. 2023, 306, 116578. [Google Scholar] [CrossRef]
  8. Herwig, T.; Wagner, W. On a robust FE2 model for delamination analysis in composite structures. Compos. Struct. 2018, 201, 597–607. [Google Scholar] [CrossRef]
  9. Terada, K.; Kikuchi, N. A class of general algorithms for multiscale analyses of heterogeneous media. Comput. Methods Appl. Mech. Eng. 2001, 190, 5427–5464. [Google Scholar] [CrossRef]
  10. Terada, K.; Hori, M.; Kyoya, T.; Kikuchi, N. Simulation of the multi-scale convergence in computational homogenization approaches. Int. J. Solids Struct. 2000, 37, 2285–2311. [Google Scholar] [CrossRef]
  11. Jänicke, R.; Diebels, S.; Sehlhorst, H.G.; Düster, A. wo-scale modelling of micromorphic continua. Contin. Mech. Thermodyn. 2009, 21, 297–315. [Google Scholar] [CrossRef]
  12. Feyel, F.; Chaboche, J.L. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials. Comput. Methods Appl. Mech. Eng. 2000, 183, 309–330. [Google Scholar] [CrossRef]
  13. Tikarrouchine, E.; Chatzigeorgiou, G.; Chemisky, Y.; Meraghni, F. Fully coupled thermo-viscoplastic analysis of composite structures by means of multi-scale three-dimensional finite element computations. Int. J. Solids Struct. 2019, 164, 120–140. [Google Scholar] [CrossRef] [Green Version]
  14. Kaczmarczyk, Ł.; Pearce, C.J.; Bićanić, N. Scale transition and enforcement of RVE boundary conditions in second-order computational homogenization. Int. J. Numer. Methods Eng. 2008, 74, 506–522. [Google Scholar] [CrossRef] [Green Version]
  15. Kaczmarczyk, Ł.; Pearce, C.J.; Bićanić, N. Studies of microstructural size effect and higher-order deformation in second-order computational homogenization. Comput. Struct. 2010, 88, 1383–1390. [Google Scholar] [CrossRef]
  16. Bacigalupo, A.; Gambarotta, L. Non-local computational homogenization of periodic masonry. Int. J. Multiscale Comput. Eng. 2011, 9, 565–578. [Google Scholar] [CrossRef] [Green Version]
  17. Tan, V.B.C.; Raju, K.; Lee, H.P. Direct FE2 for concurrent multilevel modelling of heterogeneous structures. Comput. Methods Appl. Mech. Eng. 2020, 360, 112694. [Google Scholar] [CrossRef]
  18. Liu, K.; Meng, L.; Zhao, A.; Wang, Z.G.; Chen, L.L.; Li, P. A hybrid direct FE2 method for modeling of multiscale materials and structures with strain localization. Comput. Methods Appl. Mech. Eng. 2023, 412, 116080. [Google Scholar] [CrossRef]
  19. Chen, W.J.; Tan, V.B.C.; Zeng, X.G.; Li, P. FE2 methodology for discrete cohesive crack propagation in heterogenous materials. Eng. Fract. Mech. 2022, 269, 108537. [Google Scholar] [CrossRef]
  20. Xu, J.H.; Li, P.; Poh, L.H.; Zhang, Y.Y.; Tan, V.B.C. Direct FE2 for concurrent multilevel modelling modeling of heterogeneous thin plate structures. Comput. Methods Appl. Mech. Eng. 2022, 392, 114658. [Google Scholar] [CrossRef]
  21. Zhi, J.; Poh, L.H.; Tay, T.E.; Tan, V.B.C. Direct FE2 modeling of heterogeneous materials with a micromorphic computational homogenization framework. Comput. Methods Appl. Mech. Eng. 2022, 393, 114837. [Google Scholar] [CrossRef]
  22. Yeoh, K.M.; Poh, L.H.; Tay, T.E.; Tan, V.B.C. Multiscale computational homogenisation of shear-flexible beam elements: A Direct FE2 approach. Comput. Mech. 2022, 70, 891–910. [Google Scholar] [CrossRef]
  23. Bendsøe, M.P.; Sigmund, O. Material interpolation schemes in topology optimization. Arch. Appl. Mech. 1999, 69, 635–654. [Google Scholar] [CrossRef]
  24. Xie, Y.M.; Steven, G.P. A simple evolutionary procedure for structural optimization. Comput. Struct. 1993, 49, 885–896. [Google Scholar] [CrossRef]
  25. Wang, M.Y.; Wang, X.; Guo, D. A level set method for structural topology optimization. Comput. Methods Appl. Mech. Eng. 2003, 192, 227–246. [Google Scholar] [CrossRef]
  26. Rietz, A. Sufficiency of a finite exponent in SIMP (power law) methods. Struct. Multidiscip. Optim. 2001, 21, 159–163. [Google Scholar] [CrossRef]
  27. Stolpe, M.; Svanberg, K. An alternative interpolation scheme for minimum compliance topology optimization. Struct. Multidiscip. Optim. 2001, 22, 116–124. [Google Scholar] [CrossRef]
  28. Huang, X.D.; Xie, Y.M.; Meron, B.M.C. A New Algorithm for Bi-Directional Evolutionary Structural Optimization. JSME Int. J. Ser. C Mech. Syst. Mach. Elem. Manuf. 2006, 49, 1091–1099. [Google Scholar] [CrossRef] [Green Version]
  29. Sethian, J.A.; Wiegmann, A. Structural Boundary Design via Level Set and Immersed Interface Methods. J. Comput. Phys. 2000, 163, 489–528. [Google Scholar] [CrossRef]
  30. Xia, L.; Breitkopf, P. Concurrent topology optimization design of material and structure within FE2 nonlinear multiscale analysis framework. Comput. Methods Appl. Mech. Eng. 2014, 278, 524–542. [Google Scholar] [CrossRef]
  31. Rodrigues, H.; Guedes, J.M.; Bendsøe, M.P. Hierarchical optimization of material and structure. Struct. Multidiscip. Optim. 2002, 24, 1–10. [Google Scholar] [CrossRef]
  32. Coelho, P.G.; Fernandes, P.R.; Guede, J.M.; Rodrigues, H.C. A hierarchical model for concurrent material and topology optimisation of three-dimensional structures. Struct. Multidiscip. Optim. 2008, 35, 107–115. [Google Scholar] [CrossRef]
  33. Li, Q.; Xu, R.; Wu, Q.; Liu, S. Topology optimization design of quasi-periodic cellular structures based on erode-dilate operators. Comput. Methods Appl. Mech. Eng. 2021, 377, 113720. [Google Scholar] [CrossRef]
  34. Hu, J.; Luo, Y.; Liu, S. Two-scale concurrent topology optimization method of hierarchical structures with self-connected multiple lattice-material domains. Compos. Struct. 2021, 272, 114224. [Google Scholar] [CrossRef]
  35. Hu, J.; Luo, Y.; Liu, S. Three-scale concurrent topology optimization for the design of the hierarchical cellular structure. Struct. Multidiscip. Optim. 2022, 65, 143. [Google Scholar] [CrossRef]
  36. Jantos, D.R.; Hackl, K.; Junker, P. Topology optimization with anisotropic materials, including a filter to smooth fiber pathways. Struct. Multidiscip. Optim. 2020, 61, 2135–2154. [Google Scholar] [CrossRef]
  37. Wallin, M.; Ivarsson, N.; Amir, O.; Tortorelli, D. Consistent boundary conditions for PDE filter regularization in topology optimization. Struct. Multidiscip. Optim. 2020, 62, 1299–1311. [Google Scholar] [CrossRef]
  38. Yang, A.; Wang, S.; Luo, N.; Xiong, T.F.; Xie, X.D. Massively efficient filter for topology optimization based on the splitting of tensor product structure. Front. Mech. Eng. 2022, 17, 54. [Google Scholar] [CrossRef]
  39. Ibhadode, O.; Zhang, Z.D.; Sixt, J.; Nsiempba, K.M.; Orakwe, J.; Martinez-Marchese, A.; Ero, O.; Shahabad, S.I.; Bonakdar, A.; Toyserkani, E. Topology optimization for metal additive manufacturing: Current trends, challenges, and future outlook. Virtual Phys. Prototyp. 2023, 18, e2181192. [Google Scholar] [CrossRef]
  40. Theocaris, P.S.; Stavroulakis, G.E. Optimal material design in composites: An iterative approach based on homogenized cells. Comput. Methods Appl. Mech. Eng. 1999, 169, 31–42. [Google Scholar] [CrossRef]
  41. Xia, L.; Raghavan, B.; Breitkopf, P.; Zhang, W. Numerical material representation using proper orthogonal decomposition and diffuse approximation. Appl. Math. Comput. 2013, 224, 450–462. [Google Scholar] [CrossRef]
  42. Sigmund, O.; Petersson, J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct. Optim. 1998, 16, 68–75. [Google Scholar] [CrossRef]
  43. Sigmund, O. Materials with prescribed constitutive parameters: An inverse homogenization problem. Int. J. Solids Struct. 1994, 31, 2313–2329. [Google Scholar] [CrossRef]
Figure 1. The reconstruction scheme of optimization results based on the Direct FE2 method.
Figure 1. The reconstruction scheme of optimization results based on the Direct FE2 method.
Mathematics 11 02779 g001
Figure 2. Interpolating of all element densities.
Figure 2. Interpolating of all element densities.
Mathematics 11 02779 g002
Figure 3. Iso-position filter of element densities. (a) Searching the neighbor RVEs within a circular region of radius r f i l t e r ; (b) Weight sum of the same position elements in the neighbor RVEs.
Figure 3. Iso-position filter of element densities. (a) Searching the neighbor RVEs within a circular region of radius r f i l t e r ; (b) Weight sum of the same position elements in the neighbor RVEs.
Mathematics 11 02779 g003
Figure 4. (a) Direct FE2 model and (b) DNS model of the bridge-type structure.
Figure 4. (a) Direct FE2 model and (b) DNS model of the bridge-type structure.
Mathematics 11 02779 g004
Figure 5. The RVE of Direct FE2 model and DNS model.
Figure 5. The RVE of Direct FE2 model and DNS model.
Mathematics 11 02779 g005
Figure 6. The contour plot of densities after optimization of the bridge-type structure. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Figure 6. The contour plot of densities after optimization of the bridge-type structure. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Mathematics 11 02779 g006
Figure 7. The reconstruction result of the bridge-type structure (filter radius is 60   c m ).
Figure 7. The reconstruction result of the bridge-type structure (filter radius is 60   c m ).
Mathematics 11 02779 g007
Figure 8. Normalized strain energy of the bridge-type structure.
Figure 8. Normalized strain energy of the bridge-type structure.
Mathematics 11 02779 g008
Figure 9. Comparison of the optimization result of the (a) Direct FE2 model and the (b) DNS model for the bridge-type structure.
Figure 9. Comparison of the optimization result of the (a) Direct FE2 model and the (b) DNS model for the bridge-type structure.
Mathematics 11 02779 g009
Figure 10. Reconstruction results of the bridge-type structure with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Figure 10. Reconstruction results of the bridge-type structure with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Mathematics 11 02779 g010aMathematics 11 02779 g010b
Figure 11. Topological evolution of the bridge-type structure ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 26.
Figure 11. Topological evolution of the bridge-type structure ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 26.
Mathematics 11 02779 g011
Figure 12. (a) Direct FE2 model and (b) DNS model of the MBB beam.
Figure 12. (a) Direct FE2 model and (b) DNS model of the MBB beam.
Mathematics 11 02779 g012
Figure 13. The contour plot of densities after optimization of the MBB beam. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Figure 13. The contour plot of densities after optimization of the MBB beam. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Mathematics 11 02779 g013
Figure 14. The reconstruction result of the MBB beam (Filter radius is 60   c m ).
Figure 14. The reconstruction result of the MBB beam (Filter radius is 60   c m ).
Mathematics 11 02779 g014
Figure 15. Normalized strain energy of the MBB beam.
Figure 15. Normalized strain energy of the MBB beam.
Mathematics 11 02779 g015
Figure 16. Comparison of the optimization result of the (a) Direct FE2 model and the (b) DNS model for the MBB beam.
Figure 16. Comparison of the optimization result of the (a) Direct FE2 model and the (b) DNS model for the MBB beam.
Mathematics 11 02779 g016
Figure 17. Reconstruction results of the MBB beam with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Figure 17. Reconstruction results of the MBB beam with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Mathematics 11 02779 g017
Figure 18. Topological evolution of the MBB beam ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 27.
Figure 18. Topological evolution of the MBB beam ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 27.
Mathematics 11 02779 g018aMathematics 11 02779 g018b
Figure 19. (a) Direct FE2 model and (b) DNS model of the loaded knee model.
Figure 19. (a) Direct FE2 model and (b) DNS model of the loaded knee model.
Mathematics 11 02779 g019
Figure 20. The contour plot of densities after optimization of the loaded knee model. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Figure 20. The contour plot of densities after optimization of the loaded knee model. Only the elements with the top 50 % of the density are drawn in the whole model picture. (ad) are magnified views of the local regions of four typical macroscale elements. V f 1 , V f 2 , V f 3 , and V f 4 represent the retained volume fraction of four RVEs in the macroscale element (not drawn to scale).
Mathematics 11 02779 g020
Figure 21. The reconstruction result of the loaded knee model (filter radius is 60   c m ).
Figure 21. The reconstruction result of the loaded knee model (filter radius is 60   c m ).
Mathematics 11 02779 g021
Figure 22. Normalized strain energy of the loaded knee model.
Figure 22. Normalized strain energy of the loaded knee model.
Mathematics 11 02779 g022
Figure 23. Comparison of the optimization result of the (a) Direct FE2 model and (b) DNS model for the loaded knee model.
Figure 23. Comparison of the optimization result of the (a) Direct FE2 model and (b) DNS model for the loaded knee model.
Mathematics 11 02779 g023
Figure 24. Reconstruction results of the loaded knee model with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Figure 24. Reconstruction results of the loaded knee model with different filter radius. (a) r f i l t e r = 20   c m ; (b) r f i l t e r = 40   c m ; (c) r f i l t e r = 60   c m ; (d) r f i l t e r = 80   c m .
Mathematics 11 02779 g024
Figure 25. Topological evolution of the loaded knee model ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 26.
Figure 25. Topological evolution of the loaded knee model ( r f i l t e r = 60   c m ). (a) Substep 1; (b) Substep 5; (c) Substep 10; (d) Substep 15; (e) Substep 20; (f) Substep 26.
Mathematics 11 02779 g025
Table 1. Comparison of topology optimization methods at different scales.
Table 1. Comparison of topology optimization methods at different scales.
Computational ModelScaleNumber of ElementsComputational EfficiencyNumerical Implementation Method
DNSSingle scalelargeSingle layerSoftware
FE2MultiscalesmallTwo-layer nestingCoding
Direct FE2MultiscalesmallSingle layerSoftware
Table 2. Comparison of the Direct FE2 and DNS model for the bridge-type structure.
Table 2. Comparison of the Direct FE2 and DNS model for the bridge-type structure.
Computational ModelNumber of ElementsNumber of IterationsOptimization Time
Direct FE2 307,200 26 49   m
DNS 7,680,000 46 19   h   13   m
Table 3. Comparison of the Direct FE2 and the DNS model for the MBB beam.
Table 3. Comparison of the Direct FE2 and the DNS model for the MBB beam.
Computational ModelNumber of ElementsNumber of IterationsOptimization Time
Direct FE2 307,200 27 41   m
DNS 7,680,000 40 16   h   30   m
Table 4. Comparison of the Direct FE2 and the DNS model for the loaded knee model.
Table 4. Comparison of the Direct FE2 and the DNS model for the loaded knee model.
Computational ModelNumber of ElementsNumber of IterationsOptimization Time
Direct FE2 460,800 26 50   m
DNS 11,520,000 47 64   h
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhao, A.; Tan, V.B.C.; Li, P.; Liu, K.; Hu, Z. A Reconstruction Approach for Concurrent Multiscale Topology Optimization Based on Direct FE2 Method. Mathematics 2023, 11, 2779. https://doi.org/10.3390/math11122779

AMA Style

Zhao A, Tan VBC, Li P, Liu K, Hu Z. A Reconstruction Approach for Concurrent Multiscale Topology Optimization Based on Direct FE2 Method. Mathematics. 2023; 11(12):2779. https://doi.org/10.3390/math11122779

Chicago/Turabian Style

Zhao, Ang, Vincent Beng Chye Tan, Pei Li, Kui Liu, and Zhendong Hu. 2023. "A Reconstruction Approach for Concurrent Multiscale Topology Optimization Based on Direct FE2 Method" Mathematics 11, no. 12: 2779. https://doi.org/10.3390/math11122779

APA Style

Zhao, A., Tan, V. B. C., Li, P., Liu, K., & Hu, Z. (2023). A Reconstruction Approach for Concurrent Multiscale Topology Optimization Based on Direct FE2 Method. Mathematics, 11(12), 2779. https://doi.org/10.3390/math11122779

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop