Topology Optimization for Continua Considering Global Displacement Constraint

Finding the best distribution of available material in the predetermined design domain satisfying various conditions is the target of topology optimization for continuum structures. In most topology optimization methods, the optimized goal is to find the structures with maximum stiffness. Stiffness is often closely related to the global displacement and, especially, to the maximum displacement of the structure. So a new topology method for minimizing the volume of the structure subject to global displacement is developed and a new approach to controlling the maximum displacement of the structure is proposed.


INTRODUCTION
Finding the best distribution of available material in the predetermined design domain satisfying various conditions is the target of topology optimization for continuum structures.In most topology optimization methods, the optimized goal is to find the structures with maximum stiffness.Stiffness is often closely related to the global displacement and, especially, to the maximum displacement of the structure.So a new topology method for minimizing the volume of the structure subject to global displacement is developed and a new approach to controlling the maximum displacement of the structure is proposed.
The topology optimization method based on elements is treated by dividing the design domain into finite elements and each element is taken as a design variable.The solid isotropic material and penalization (SIMP) method [1] and [2] transfers design variables from the discrete variables to contiguous ones with penalization.The evolutionary structural optimization (ESO) method [3] and [4], and its later version, the bi-directional ESO (BESO) method [5], remove inefficient material from the structure based on certain predefined criteria.The level set method represents the structure using a level set model which is embedded in a scalar function.Rong and Liang [6] and Wang et al. [7] investigated the level set method for topology optimization.However, the element-wise topology optimizations exhibit various numerical problems, such as grey-scale, checkerboard pattern and mesh dependency.Therefore, topology optimization methods based on nodal design variables were developed to avoid these problems.Matsui and Terada [8] proposed the concept of continuous distribution.A Q4/Q4 continuum structural topology optimization is investigated by Rahmatalla and Swan [9].
A 3-D structural topology optimization and novel surface-smoothing scheme based on SIMP and subelement bilinear interpolation was developed using node densities by Song and Kim [10].An adaptive density point refinement approach for continuum topology optimization on the basis of an analysismesh separated material density field description based on nodal design variables was presented by Wang et al. [11].Wang et al. [12] proposed topological optimization of structures using a multilevel nodal density-based approximant.The nodal density field using the non-local Shepard function method is transformed to a practical elemental density field via a local interpolation with the elemental shape function.
The presented topology optimization method is based on nodal densities and utilizes the rational approximation for material properties (RAMP) interpolation scheme proposed by Stolpe and Svanberg [13].The discrete nodal topological variables ρ i that only take value 0 or 1 are replaced by continuous topological variables between 0 and 1.Consequently, the difficulty of the discrete optimization is avoided by penalization.It is assumed: where ρ x is the density of point x, i.e. ρ(x), v is the penalization factor (v is a parameter which in some sense corresponds to p in the SIMP approach), which makes the intermediate densities approach either 0 (void) or 1 (solid).Here, v = 5 is set.The relation between the Young's modulus and the material density at point x is expressed by: where E 0 is Young's modulus of the full solid material.The function f(ρ x ) has the following properties: The volume of an element is given by: where V i is the volume of the i th element, V i 0 is the original volume of the i th element.
In this study, minimum volume with a reference domain Ω in R 3 is considered while satisfying the global displacement constraint for the structure.
where V is the structural volume being optimized, u j f is the displacement of the j th degree of freedom of the structure under the f th load case, N dof is the total degrees of freedom of the structure, U f is its constraint limit, L is the number of the load cases acting on the structure, and ρ i is the density of the i th node, ρ i is its lower limit, N n is the total number of nodes.Here, the small positive lower bound ρ i = 0.0001 is set so that the structure optimized is always kept non-singular in the optimization process.

EQUIVALENT MAXIMUM DISPLACEMENT
When the maximum displacement of the structure does not exceed a specified value, the global displacement constraint is satisfied.So, the maximum displacement is naturally the ideal design criterion of optimization models.However, the location of the maximum displacement usually varies following the change of material distribution in the optimization process, so the maximum function is not differentiable.To solve this problem,the maximum displacement needs to be smoothed, and the p-norm or the Kreisselmeier-Steinhauser (KS) function could be used.Similar to the p-norm, geometric average displacement (GAD) is introduced by Kreisselmeier [14] and Qiao and Liu [15], which is expressed as: Here, N e is the total number of elements, Ω denotes the volume of the design region, p is the norm parameter and Ω e is the element e solid volume and the displacement of any point inside the element e is δ.
On one extreme, as the displacement norm parameter p approaches infinity, the U a approaches the maximum displacement, and there is no added smoothness.On the other extreme, when p approaches 1, there is excessive smoothness but U a approaches the average displacement.A good choice for p should, therefore, provides adequate smoothness so that the optimization algorithm performs well and an adequate approximation of the maximum displacement value so that the optimized design satisfies the imposed displacement constrains.
Only when p approaches infinity, can the U a arrive at the constraint of the maximum displacement.To remedy this deficiency, a normalized global displacement measure is proposed to better approximate the maximum displacement.The normalized p-norm displacement uses information from the previous optimization iteration to scale, and the p-norm displacement as α|U a |, so that it better approximates the maximum displacement.The maximum displacement max(|u| k-1 ) and the p-norm displacement values from the previous optimization iteration k -1 are used to define our evolving normalized p-norm displacement constraint at each iteration k as: where α is calculated at each optimization iteration k ≥ 1 as follows: . Note that the constraint α|U a | ≤ U f is nondifferentiable because the value of α is changing in a discontinuous manner and results in a slightly different optimization problem at every iteration.However, as the optimization converges the changes between successive design iterations diminish and hence α converges, thereby reducing the effects of the nondifferentiability and inconsistency.

NODAL DESIGN VARIABLES AND THE INTERPOLATION SCHEME
Element-independent nodal densities are the design variables in this method.Thus, the relative density at any point is interpolated with an interpolation scheme, which determines the topology, stiffness and volume of material.

Identifying the Nodes
As proposed by Guest et al. [16] and Kang and Wang [17], the scale parameter r min is set to identify the nodes that influence the density of point x.Nodes are included in the influence domain if they are located within a distance r min of the point x.This can be visualized by drawing a sphere of radius r min centered at the point x, thus generating the spherical subdomain Ω x .Nodes located inside Ω x contribute to the computation of density ρ(x) of point x.As the mesh is refined, r min and consequently Ω x do not change.The only difference between the two meshes is the number of nodes located inside Ω x , and included in the interpolation function.This is essential to generate mesh-independent solutions.

Shepard Interpolation Scheme
Interpolation provides a continuum of density field and mesh-independence, which might alleviate numerical instability and checkerboard effects [18].In implementing continuum structural topology optimization formulations, many functions are available to interpolate nodal density onto the points inside the element space; -for example, the standard C 0 shape functions used in the finite element method.However, each node's shape function influences only the elements connected to that node.Mesh independency cannot be obtained when interpolation functions are influenced by mesh size.They should be based on a physical length scale that does not change following mesh refinement.
Shepard interpolation is proposed by Shepard [19], and used by Kang and Wang [17] to achieve mesh independency.Let ρ i (i = 1, 2, ..., n) denote a set of density of nodes inside Ω x at the associated point x = (X, Y, Z), where (X, Y, Z) define the point x location in the Cartesian coordinate system.Thus, the relative density at point x is interpolated by the nodal densities inside the influence domain Ω x with Shepard's interpolation method.
where S x is the sub-domain of design variable located within the influence domain Ω x of point x, and ρ i is the density value of the i th node.x i is the position of the point associated with the i th node.The corresponding interpolation function ND i (x) is defined as: where R i (x) = 1 / [r i (x)] and r i (x) = |x -x i | 2 being the Euclidean distance from the points x to x i .n is the number of nodes inside the influence domain Ω x .
In the method of element independent nodal variable density, the density in the element space is not constant, and the global density field of the structure has C 0 continuity.It is easy simple to know from the bounded property of the Shepard interpolation that 0 ≤ ρ x ≤ 1 holds if 0 ≤ ρ i ≤ 1 (i ∈ S x ).Moreover, the property ND i (x) ≥ 0 also guarantees that the derivative of the density with respect to the design variable will be always non-negative.This property is essential to guaranteeing a correct searching direction in seeking the optimal material distribution by a gradient-based algorithm.

SENSITIVITY ANALYSIS
The solution of the gradient-based optimization problem requires the computation of sensitivities of the objective function and the constraints.In a finite element analysis, the static behavior of a structure for any load case can be expressed by the following equilibrium equation: where K is the global stiffness matrix of a structure being optimized and, U and F are the global nodal displacement and nodal load vector, respectively.
In the finite element method, the displacement of any point δ can be expressed by: Here, N, u denote the general shape function matrix and displacement vector of the i th element, respectively.
The adjoint method can be used to determine the sensitivity of displacement by introducing a vector of Lagrange multiplier λ.The modified p-norm displacement can be expressed as: where the term λ T (KU -F) is equal to zero and, therefore, the modified displacement is identical to the original one.Taking derivatives of Eg. ( 12) with respect to the design variable ρ gives: The Eq. ( 13) can be rewritten as: 1 To eliminate the unknown ∂U / ∂ρ from the sensitivity expression, let: So the adjoint vector is defined as: Here, P λ denotes the adjoint load of adjoint Eq. for yielding the adjoint vector λ.
Thus, the sensitivity of the p-norm displacement function is: where N c is the number of element influenced by ρ i .K 0 e is the element stiffness matrix of the e th element of the solid material.
The sensitivity of the p-norm displacement requires the computation of the sensitivity of the stiffness matrix with respect to the design variable.The derivative of the elemental stiffness matrix with respect to the design variable is expressed by: where B is the usual displacement-strain matrix and D 0 corresponds to the constitutive matrix of the solid material.
For example, the formulation of the constitutive matrix for 3D isotropic solid structures is: Since the stiffness matrix integrand is evaluated at the Gauss points, the densities at these Gauss points are directly computed from the design variables using the interpolation function.Numerical quadrature, such as Gaussian quadrature, is commonly reduced to the evaluation and summation of the stiffness integrand at specific Gauss points.
The sensitivity analysis of the objective function in Eq. ( 4) is calculated similar to that of Eq. ( 19).The derivative of the total material volume with respect to the design variables can be computed by the Gaussian quadrature method over the influence-domain.
So the total volume is: The derivative of V with respect to the design variable:

TOPOLOGICAL OPTIMIZATION WITH VARYING DISPLACEMENT LIMITS
In order to make the approximation functions of the displacement constraint in Eq. ( 4) hold true at each iteration, an equivalent optimization model (Eq.( 25)) with varying displacement constraint limits is built.
At each iteration, a quadratic approximation model of the true objective function that satisfies the Taylor expansion is built around the current point.The model is assumed to be a good representative of the objective function in a so-called trust region [20].Trust regions are used to ensure the robustness of the iteration and make progress toward feasibility and optimality [21].
where U l f is expressed as: In the above, β is a displacement limit changing factor.Typical values of β between 0.01 and 0.20 have been used for displacement constraints in the example problems in this paper.U a f k is the p-norm displacement of the structure under the f th load case.
are varied by Eq. ( 26) at every iteration.
Assuming that only t i = 1/ρ i is changed and is treated as a variable, the first-order series expansion for the p-norm displacement U a at t i (i = 1, 2, ..., N nod ) can be expressed as: where N n is the number of the nodes in the design domain.
Thus, the p-norm displacement in the next iteration, U a f k+1 , can be estimated by the current iteration k.
(25) can be transferred into Eq.( 28): min : . .( ) If the constant items in the objective function are omitted, solving Eq. ( 28) can be transferred to solving Eq. ( 29): The programming solving the problem of Eq. ( 29) can be transferred into solving dual programming problem by using the dual theory [22] and [23].

NUMERICAL EXAMPLES
This section illustrates the proposed approach with 3D applications.For simplicity, all the quantities are dimensionless.In addition, Young's modulus is chosen as 2.1×10 11 and Poisson's ratio as 0.3 for all examples.Let d denote the length of cuboid element diagonal.In the following examples, r min is set to 1.5d while keeping the same displacement mesh size.As the computing platform, we have used a personal computer with the commercial software package ABAQUS has been used for FEA in this study.
Fig. 1a shows a 3D cantilever beam with length of 8, height of 5, and width of 2. The beam is fixed at the left end.A concentrated load P = 4800 is applied downward at the center of the right end.The initial displacement at the constraint point, the center of the right end, is 2.851×10 -6 .A displacement constraint is taken as 7.0×10 -6 .The cantilever domain is discretized into a mesh of 40×25×10 B8 elements which results in a total of 10,000 elements.There are 41×26×11 design variable points distributed within the initial design domain.The optimized topology is shown in Fig. 1b.Then, the same problems are solved using the proposed method with varying p parameter in which p is set to even numbers (4, 10, 20, 40) as shown in Fig. 2. With the increase of p parameter, the displacement constraint converges to the same condition (Fig. 3), the material distribution (Fig. 2) and the volume (Fig. 4) move close to the convergence.Fig. 5 shows a 3D, simply supported beam with a length of 4, height of 1, and width of 0.4.An external force P = 4000 is applied to the center of top area.The initial displacement at the constraint point, the center of the top area, is 1.1×10 -7 .A displacement constraint is taken as 2.0×10 -7 .The design domain is discretized into mesh size of 80×20×8 B8 elements.p is set to 20.When the nodal density values are used to determine a smooth iso-line that describes the boundary of the optimization layout, as a result, a smooth optimal topology can be obtained in Fig. 6a.Its final volume is 0.947.The results obtained from the element-based approach is shown in Fig. 6b.Its final volume is 1.062.When the mesh is coarser, and the mesh size of 50×12×5 is used, the results obtained from the element-based approach and the proposed approach based on element independent nodal density are shown in Figs.7a and 3b, respectively.Their final volumes are 1.071 and 1.155.
These figures show that, for the different displacement mesh size, the topology obtained from element independent nodal density has a much better resolution and is smoother than that of the elementbased approach.The solution attained by the proposed method exhibits no checkerboard patterns or mesh dependency.
When the same mesh is used, the computational cost for the topology optimization based on element independent nodal density is higher than the elementbased approach.This is mainly attributable to the large number of density nodes in the influence domain.However, the topology resolution resulting from the proposed approach based on the proposed method is higher than that of the element-based approach.To improve the efficiency of the proposed approach, especially for a 3D large-scale optimization problem, the parallel programming technique could be used to carry out the finite element analysis and the optimization procedure.
The total CPU time and the CPU time spent on the sensitivity analysis in every optimization iteration are 144 and 123 seconds, respectively.When the codes are reprogrammed with the OpenMP and four threads are used, the total CPU time and the CPU time spent on the sensitivity analysis decrease to 50 and 42 seconds, respectively.By taking this approach, it is possible to obtain benefits from parallelization without the need for extensive modification to the code structure.

CONCLUSIONS
This paper has developed a topology optimization method for minimizing the volume of a structure subject to the global displacement constraint.In contrast to the element-based procedure, here we take the nodal density as the design variable, which is interpolated into any point by Shepard functions.This technique avoids checkerboard patterns and meshdependency for low order finite elements.With the help of the global displacement constraint, an optimal structure with appointed deformation can be obtained, and it is unnecessary to know where the maximum displacement is.The proposed method is highly useful with regard to practical engineering applications.The numerical examples demonstrate the effectiveness of the proposed method with respect to the optimal solution and convergence.

Fig. 1 .
The 3D cantilever beam; a) design domain and boundary conditions; b) topology optimization result obtained from the proposed method