Homogenization Topology Optimization Method Based on Continuous Field

Abstract

In order to overcome numerical instabilities such as checkerboards, meshdependence in topology optimization of continuum structures, a new implementation combined with homogenization method without introducing any additional constraint parameter is presented. To overcome the shortcoming of continuous material distribution by the introduction of finite element approximation, moving least square or modified filter functions are adopted as interpolation function. The method can be viewed as a nature extension of node-based homogenization method and named as material point homogenization method. Continuous size field and continuous density field are constructed, and structural responses' sensitivities are derived. Several representative numerical examples are presented to demonstrate the capability and the efficiency of the proposed approach against some classes of numerical instabilities.

1. Introduction

Topology optimization refers to optimal design problems in which the topology of the structure is allowed to change in order to improve the performance of the structure. Relative to size optimization and shape optimization, this optimization class is regarded as one of the most challenging optimization problems in the field of structural optimization. Much research has been devoted to topology optimization over the last decades. Starting with the landmark paper of Bendsøe and Kikuchi [1], numerical methods for topology optimization have been investigated extensively since the late 1980s. The area of topology optimization of continuum structures is now dominated by methods that employ the material distribution concept. The typical ones are the homogenization and variable density method. In the homogenization method, the material property of each design cell is computed by the homogenization theory, and the optimal topology is obtained by solving a material distribution problem. However, the homogenization method may not yield the intended results with infinitesimal pores in the materials that make the structure not manufacturable.

An engineering alternative approach called the variable density method was introduced in [2]. Instead of using the homogenization theory to evaluate the material property for design cells, this approach assumes some explicit relationships between the density and the material property without considering their microstructures. This approach is very attractive to the engineering community because of its simplicity. The power-law approach is physically permissible as long as simple conditions on the power are satisfied (e.g., $p\ge 3$ for Poisson's ration equal to 1/3). Nowadays, this method has been widely developed to problems with multiple constraints, multiple physics, and multiple materials.

A simple method for shape and layout optimization, called “evolutionary structural optimization” (ESO), has been proposed by Xie and Steven which is based on the idea of gradually removing inefficient material to achieve an optimal design [3]. Following this basic approach, there have been a number of modifications and refinements such as BESO, where not only are elements removed but are also added in high stressed areas [4]. This kind of method was developed for various problems of structural optimization including stress considerations, stiffness constraints, frequency optimization, and buckling problems. Nevertheless, ESO has some drawbacks and weaknesses. The first is its expensive computational cost and the second is questionable validity proposed by Zhou and Rozvany [5, 6].

Most material distributed methods suffer from the checkerboard patterns, mesh dependencies, and grey-scale problem. As reviewed by Sigmund and Petersson [7], a large number of works have suggested different regularization schemes to be used in connection with topology optimization. These methods can roughly be divided into fours categories: (1) constraint methods such as perimeter control, global or local gradient control [8–10]; (2) the use of higher-order elements, nonconforming or hybrid elements; (3) various mesh-independent filtering [11]; (4) other methods like wavelet space [12]. These schemes strongly rely on the artificial parameters for additional constraints to avoid the instabilities, and there are no rational guidelines to a priori determine the best parameter. In all element-based computational framework, the initial design space is always discretized by uniform finite elements, and the design variables are assumed to be constant within each finite element. Although element-based computational framework is efficient in computation and has been applied successfully for the solutions of many industrial optimization problems, it still has some undesirable features. In element-based model, the finite element grid is used both for representing the topology of the structure and performing physical analysis. Matsui and Terada [13] first proposed a checkerboard-free homogenization method by virtue of the continuous FE approximation of design variables. Analogously, Ramatalla and Swan [14] developed a so-called nodal variable density method which ensures the continuity of density field. Anther modified Q4/Q4 element is presented to demonstrate the efficiency in generating topologies with high resolution [15]. Belytschko et al. [16] developed level set method in terms of the nodal variables that control an implicit function description of the shape. All above methods based on nodal design variables can overcome checkerboard patterns without any restriction methods. Very recently, meshfree method combined with nodal or Gauss point density as design variable is explored to cope with topology optimization problems with linear elasticity or geometric nonlinearity [17].

In this paper, we propose a checkerboard-free and mesh-independence topology optimization method without introducing any additional constraint scheme. This aim is accomplished by the introduction of continuous material distribution. By virtue of this continuous approximation of design, discontinuous distribution like checkerboard pattern disappears without any schemes. Representative numerical examples are presented to demonstrate the capability of the proposed method.

2. Homogenization Theory for Media with Periodic

The computation of effective properties plays a key role for the topology optimization [18–20]. Suppose that a periodic microstructure is assumed in the neighborhood of an arbitrary point x of a given linearly elastic structure. The periodicity is represented by a parameter δ which is very small, and the elastic tensor $E_{ijkl}^{\delta }$ is given in the form [21]

$E_{ijkl}^{\delta }=E_{ijkl}\left(x,\frac{x}{\delta }\right),$

(1)

where $y\to E_{ijkl}(x,y)$ is Y-periodic, with cell $Y=[Y_{1R},Y_{1L}]\times [Y_{2R},Y_{2L}]$ of periodicity. Here, x is the macroscopic variation of material parameters, while x/δ gives the microscopic, periodic variations. Now, suppose that the structure is subjected to a macroscopic body force and a macroscopic surface traction. The resulting displacement field $u^{\delta }(x)$ can then be expanded as

$u^{\delta }\left(x\right)=u_0\left(x\right)+\delta u_1\left(x,\frac{x}{\delta }\right)+\dots ,$

(2)

where the leading term $u_0(x)$ is a macroscopic deformation field that is independent of the microscopic variable y. It turns out that this effective displacement field is the macroscopic deformation field that arises due to the applied forces when the rigidity of the structure is given by the effective rigidity tensor

$E_{ijkl}^H\left(x\right)=\frac{1}{\left|Y\right|}{{\displaystyle \int }}_Y\left[E_{ijkl}\left(x,y\right)-E_{ijpq}\left(x,y\right)\frac{\partial {\chi }_p^{kl}}{\partial y_q}\right]dy.$

(3)

Here, ${\chi }^{kl}$ is a microscopic displacement field that is given as the Y-periodic solution of the cellproblem in weak form

${{\displaystyle \int }}_Y{E}_{ijpq}\frac{\partial {\chi }_p^{kl}}{\partial y_q}\frac{\partial v_i\left(y\right)}{\partial y_i}dy\hspace{0.17em}=\hspace{0.17em}\frac{1}{\left|Y\right|}{{\displaystyle \int }}_Y{E}_{ijpq}\frac{\partial v_i\left(y\right)}{\partial y_j}dy,\hspace{1em}\forall v∊U_Y,$

(4)

where U _Y denotes the set of all Y-periodic virtual displacement fields. From (3), we see that the effective module can be computed by solving three analysis problems for the unit cell.

In order to carry out topology optimization, numerical method solving the homogenization equations is adopted. Hexagonal microstructure and finite element model for honeycomb base cell are shown in Figures 1 and 2. This microstructure is nearly isotropic mechanical behavior in macroscale. Here, we have assumed that the microstructure is composed of the material that has Young's modulus E = 210 GPa and Poisson's ration ν = 0.3. Regularized density can be formulated as

${\rho }^H={\rho }_0\left[1-{\left(1-r\right)}^2\right].$

(5)

Figure 1

Hexagonal microstructure for topology optimization.

Figure 2

Finite element model of honeycomb base cell.

To find the elements of the homogenized elasticity matrix in function form from the discrete values obtained at the sampling points, least squares approximation polynomials were employed here. Several of the elements of the elasticity matrix with respect to the dimensions of the holes are illustrated in Figure 3.

Figure 3

Material properties as functions of regularized density.

Figure 4

Illustration of the support domain.

3. Construction of Continuous Field

3.1. Construction of Continuous Field

Topology optimization methods based on continuous distribution of material densities or equivalent design variables can overcome checkerboard patterns by virtue of the continuous approximation. However, some shortcomings are pointed out in the literature [22]. First, the material field is continuous but not smooth based on low-order element shape function interpolation. Second, structural analysis method should be finite element method, and the number and location of nodes are limited by this discrete grid. The use of moving least square approximation (MLS) instead of finite element shape function interpolation is adopted to overcome these difficulties above. To succeed to this idea, this paper assume that the microstructures are continuously distributed everywhere in the fixed design domain. The design variables or densities are continuous functions of position representing the material point located in the design domain. Here, the honeycomb cell size is used as design variables. Two kinds of continuous fields are presented in this paper. These are continuous size distribution field and continuous density distribution field. In continuous size distribution field, honeycomb cell size at any location $r^e$ in the continuum structure is interpolated as

$r^e={{\displaystyle \sum }}_{i=1}^m{\Phi }_i{r}_i,$

(6)

where m is the total number of topological variables in domain of influence and Φ _i the interpolation shape function. r _i is size design variable. Then the density at this point is the funcion of r _e and can be derived as follows:

${\rho }^e=1-{\left(1-r^e\right)}^2.$

(7)

In continuous density distribution field, the density for design variable r _i can be given by

${\rho }_i=1-{\left(1-r_i\right)}^2.$

(8)

So density at any location in the continuum structure is interpolated as

${\rho }^e={{\displaystyle \sum }}_{i=1}^m{\Phi }_i{\rho }_i.$

(9)

In (6) and (9), if Φ _i is taken as linear or bilinear shape function for interpolation in two-dimensional problems, the method is the node-based style method. To overcome shortcomings in node-based style method, three different shape functions are taken as interpolated functions in this paper. Not only MLS approximation is adopted, but also two additional shape function-based filtering schemes are listed as follows. Take the continuous size distribution field as an example

$\begin{array}{c}{\Phi }_i=\frac{\omega \left(r^j\right)}{{\sum }_{i∊{\Omega }^j}\omega \left(r^j\right)},\end{array}$

(10)

$\begin{array}{c}\omega \left(r^j\right)=\frac{\left[r_{\min \hspace{0.17em}}-\mathrm{dist}\hspace{0.17em}\left(i,j\right)\right]}{r_{\min \hspace{0.17em}}}\hspace{1em}\text{or}\hspace{1em}\omega \left(r^j\right)=\frac{\mathrm{dist}\hspace{0.17em}\left(i,j\right)}{r_{\min \hspace{0.17em}}},\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\left\{i∊N|\mathrm{dist}\hspace{0.17em}\left(i,j\right)r_{\min \hspace{0.17em}}\right\},\end{array}$

(11)

where r _min is the filtering radius and $\mathrm{dist}(i,j)$ is the distance of material point i and j. This expression is widely used to eliminate numerical instabilities in topology optimization. It should be emphasized that when MLS or (10) are adopted as interpolation scheme, the location and number of design variables are not limited by discrete grid. If the method in literature [13] is named as node-based homogenization method, the method proposed in this paper can be named as material point homogenization method. According to (11), the filtering radius can control the minimum size of structural members and minimum size of holes. This paper named these two different manners as minimum size and minimum hole interpolation, respectively [23, 24].

3.2. Structural Response and Sensitivity Analysis in Continuous Field

For convenient purpose, finite element method is adopted as structural analysis tool in this present study. Global stiffness matrix and total volume are expressed as

$\begin{array}{c}\mathbf{K}={\displaystyle \sum }_{e=1}^M{\mathbf{k}}_e={\displaystyle \sum }_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}{\mathbf{B}}^{\text{T}}\mathbf{D}\left({\rho }^e\right)\mathbf{B}\hspace{0.17em}d\Omega ,\\ V={\displaystyle \sum }_{e=1}^M{V}_e={\displaystyle \sum }_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}\rho \left({\rho }^e\right)d\Omega .\end{array}$

(12)

Here M is the total number of finite elements. B is the strain matrix and D is the elastic matrix. In numerical integral, ρ ^e is evaluated according to (7) or (9) at Gauss point.

The formulations for structural sensitivities are a little different in two different continuous fields. For continuous size distribution field, the derivatives of the stiffness matrix and volume with respect to variable r _i can be written as

$\begin{array}{c}\frac{\partial \mathbf{K}}{\partial r_i}=\frac{\partial \mathbf{K}}{\partial {\rho }^e}·\frac{\partial {\rho }^e}{\partial r^e}·\frac{\partial r^e}{\partial r_i}={{\displaystyle \sum }}_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}2{\Phi }_i\left(1-r^e\right){\mathbf{B}}^{\text{T}}\frac{\partial \mathbf{D}\left({\rho }^e\right)}{\partial {\rho }^e}\mathbf{B}\hspace{0.17em}d\Omega ,\\ \frac{\partial V}{\partial r_i}=\frac{\partial V}{\partial {\rho }^e}·\frac{\partial {\rho }^e}{\partial r^e}·\frac{\partial r^e}{\partial r_i}={{\displaystyle \sum }}_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}2{\Phi }_i\left(1-r_e\right)\frac{\partial \rho \left({\rho }^e\right)}{\partial {\rho }^e}d\Omega .\end{array}$

(13)

For continuous density distribution field, the derivatives of the stiffness matrix and volume with respect to variable r _i can be written as

$\begin{array}{cc}\hfill & \frac{\partial \mathbf{K}}{\partial r_i}=\frac{\partial \mathbf{K}}{\partial {\rho }^e}·\frac{\partial {\rho }^e}{\partial {\rho }_i}·\frac{\partial {\rho }_i}{\partial r_i}\hfill \\ \hfill & \hspace{1em}=2\left(1-r_i\right)\left({{\displaystyle \sum }}_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}{\Phi }_i{\mathbf{B}}^{\text{T}}\frac{\partial \mathbf{D}\left({\rho }^e\right)}{\partial {\rho }^e}\mathbf{B}\hspace{0.17em}d\Omega \right),\hfill \\ \hfill & \frac{\partial V}{\partial r_i}=\frac{\partial V}{\partial {\rho }^e}·\frac{\partial {\rho }^e}{\partial {\rho }_i}·\frac{\partial {\rho }_i}{\partial r_i}=2\left(1-r_i\right)\left({{\displaystyle \sum }}_{e=1}^M{{\displaystyle \int }}_{{\Omega }_e}{\Phi }_i\frac{\partial \rho \left({\rho }^e\right)}{\partial {\rho }^e}d\Omega \right).\hfill \end{array}$

(14)

4. Construction and Solution for Topological Model

A structure with maximum global stiffness provides a minimum value for the external work with the real displacement field or minimum mean compliance. Therefore, the structural optimization problem with the mean compliance as the objective function can be constructed as

$\begin{array}{cc}\hfill \text{min:}& \hspace{1em}C={\mathbf{F}}^{\text{T}}\mathbf{U}\hfill \\ \hfill \text{s}.\text{t}.& \hspace{1em}\frac{V}{V_0-f}0,\hfill \\ \hfill & \hspace{1em}\mathbf{K}\mathbf{U}=\mathbf{F},\hfill \\ \hfill & \hspace{1em}{\underset{\_}{r}}_ir_i\overline{r_i}\hspace{1em}\left(i=\mathrm{1,2},\dots ,\text{N}\right).\hfill \end{array}$

(15)

Here, C is the total compliance. U and F are the global displacement and force vectors, respectively. K is the global stiffness matrix. V and V₀ are the material volume and design domain volume, respectively. ${\underset{\_}{r}}_i$ is the minimum size to avoid singularity. ${\overline{r}}_i$ is the maximum size value. To avoid the singularity when design variable takes the value r = 1.0, a modified heuristic updating scheme for the design variables if a single volume constraint is adopted here [23]

$r_i^{\left(k+1\right)}=\{\begin{array}{c}\max \hspace{0.17em}\left(r_{\min \hspace{0.17em}},\left(1-m\right)r_i^{\left(k\right)}\right),\hfill \\ {\left(B_i^{\left(k\right)}\right)}^{\beta }{r}_i^{\left(k\right)}\max \hspace{0.17em}\left(r_{\min \hspace{0.17em}},\left(1-m\right)r_i^{\left(k\right)}\right)\hfill \\ {\left(B_i^{\left(k\right)}\right)}^{\beta }{r}_i^{\left(k\right)},\hfill \\ \max \hspace{0.17em}\left(r_{\min \hspace{0.17em}},\left(1-m\right)r_i^{\left(k\right)}\right)\hfill \\ \hspace{1em}<{\left(B_i^{\left(k\right)}\right)}^{\beta }{r}_i^{\left(k\right)}<\min \hspace{0.17em}\left(1,\left(1+m\right)r_i^{\left(k\right)}\right)\hfill \\ \min \hspace{0.17em}\left(0.999,\left(1+m\right)r_i^{\left(k\right)}\right),\hfill \\ \min \hspace{0.17em}\left(0.999,\left(1+m\right)r_i^{\left(k\right)}\right){\left(B_i^{\left(k\right)}\right)}^{\beta }{r}_i^{\left(k\right)},\hfill \end{array}$

(16)

where superscript k is the number at the kth iteration, move limit m and numerical damping coefficient η are adopted to ensure robustness of the algorithm. B _i is found from the optimality condition as

$B_i=\frac{\max \hspace{0.17em}\left(0,-\partial C/\partial r_i\right)}{l\left(\partial V/\partial r_i\right)},$

(17)

where l is the Lagrangian multiplier that can be found by a bisectioning algorithm. The sensitivities of objective function $\partial C/\partial t_i$ and volume function $\partial V/\partial t_i$ can be computed according to (13)–(17). The process of optimization stops if the following convergence condition is satisfied:

$\frac{\Vert {\mathbf{r}}^{\left(\kappa +1\right)}-{\mathbf{r}}^{\left(\kappa \right)}\Vert }{\Vert {\mathbf{r}}^{\left(\kappa +1\right)}\Vert }1\text{\%}.$

(18)

5. Numerical Examples

It is well known that the use of higher-order finite elements is effective against checkerboard patterns and for clearer topologies though much computational costs are required. To illustrate the effect of eliminating the numerical instabilities, bilinear quadrangle and linear triangle plane stress elements are adopted for discretization. Different kinds of plane stress elements and their nodes and Gauss points are illustrated in Figure 5. In every numerical example, 50% volume constraint with compliance minimization is adopted. The implementation is done in MATLAB2007.

Illustration for plane element.

Example 1.

It is the optimal topological design of the so-called MBB-beam, which has been used to be benchmark in many literatures. The design domain and its dimensions are shown in Figure 6. Due to symmetry, we model only half the design domain, which is discretized with 60 by 20 bilinear quadrilaterals. Optimum topology using continuous size field or continuous density field is obtained in Figures 7 and 8.

Figure 6

Design domain and boundary condition for MBB beam.

Optimum topology using continuous size field.

Optimum topology using continuous density field.

It is obvious that final topology is with many discontinuous scattered points in low-order element shape function. This phenomenon is similar to the checkerboard patterns of using element-based topological variables. Therefore, the low-order element shape function interpolation cannot ensure the elimination of appearing discontinuous scattered points. On the contrary, less discontinuous scattered points appear in MLS approximation interpolation. And ideal optimum topological configurations are obtained using minimum size or hole interpolation.

Example 2.

The geometry loading and boundary conditions of the short cantilever beam (SCB) are illustrated in Figure 9. A point load is applied at the bottom point on the right hand side. This short cantilever beam is still used as typical example in topology optimization. The design domain is discretized with 96 by 30 triangular elements. Optimum topologies using continuous size field or continuous density field in different loading cases are obtained in Figures 10 and 11.

Figure 9

Design domain and boundary condition for SCB beam.

Optimum topology using continuous size field for case 2.

Optimum topology using continuous density field for case 2.

6. Conclusion

In this paper, we have proposed the method of continuous approximation of material distribution for topology optimization, which is consistent with the mathematical models derived from the homogenization method. The spatial interpolation by low-order element shape function cannot ensure that the optimum topologies are free from checkerboard patterns and mesh dependence. To overcome this shortcoming, MLS approximation and modified filter functions are adopted as interpolation function to construct continuous size or continuous density field. The method can be viewed as a nature extension of node-based homogenization method. Several typical two-dimension numerical examples are used to test the methods. Results are shown to demonstrate the feasibility and generality of the method proposed in this paper.