A Cell-Based Linear Smoothed Finite Element Method for Polygonal Topology Optimization

Changkye Lee,Sundararajan Natarajan,Seong-Hoon Kee and Jurng-Jae Yee,★

1University Core Research Center for Disaster-Free&Safe Ocean City Construction,Dong-A University,Busan,49315,Korea

2Department of Mechanical Engineering,Indian Institute of Technology Madras,Chennai,600036,India

3Department of Architectural Engineering,Dong-A University,Busan,49315,Korea

ABSTRACT The aim of this work is to employ a modified cell-based smoothed finite element method(S-FEM)for topology optimization with the domain discretized with arbitrary polygons.In the present work,the linear polynomial basis function is used as the weight function instead of the constant weight function used in the standard S-FEM.This improves the accuracy and yields an optimal convergence rate.The gradients are smoothed over each smoothing domain,then used to compute the stiffness matrix.Within the proposed scheme,an optimum topology procedure is conducted over the smoothing domains.Structural materials are distributed over each smoothing domain and the filtering scheme relies on the smoothing domain.Numerical tests are carried out to pursue the performance of the proposed optimization by comparing convergence,efficiency and accuracy.

KEYWORDS Smoothed finite element method;linear smoothing function;topology optimization;solid isotropic material with penalization (SIMP);polygonal finite element cell

1 Introduction

The structural topology optimization is playing the role as one of the effective tools to seek an optimal material distribution within a given domain showing the best structural performance [1].Since the inception of homogenization method by Bendsøe et al.[2],a variety of topology optimization schemes have been introduced,viz.solid isotropic material with penalization (often called SIMP) [3–5],evolutionary structural optimization (ESO) [6] and level-set method [7,8],to name a few.Among these approaches,SIMP technique,in particular,is popular due to its ease of implementation;thus,it embraces a wide area of engineering fields [9].To tackle a specific design domain in engineering fields,a reliable discretization of the domain is required,which is supplemented with design variables and appropriate constraint conditions.Due to its robustness,the finite element method (FEM) is employed to structural optimization.In FEM,typically lowerorder triangular or quadrilateral cells are used.A triangular finite element cell is easy to use but its performance is not high as quadrilaterals.On the other hand,a quadrilateral cell provides relatively better results compared to triangles;yet,it often meets difficulty to generate meshes for very complex problem domains.In addition,it is known that quadrilateral finite elements may cause a rotational mesh-dependency [10].To alleviate these difficulties,one viable alternative is to relax the mesh topology constraint,by employing element with arbitrary number of sides within the framework of the polygonal finite element method.

Not only polytopes show great flexibility in meshing the geometrically complex-shaped problem domains,but also are able to produce more accurate solutions than the standard triangular/quadrilateral finite element cells [11].An attempt to topology optimization with polygonal cells is introduced by Talischi et al.[10].Then,they developed a Matlab code of the polygonal mesh generator for a general purpose [12] to deliver an effective polygonal topology optimization [13,14].Polygonal topology optimization is extended to fluid flow optimization [15],dynamics [16] and polytree-enhancement [11].

Meanwhile,an endeavour to improve the quality of triangular cells has undertaken to bring its convenience of mesh generation.Strain smoothing approach (so called S-FEM) [17–19] has received notable attention as a reliable candidate to overcome difficulties of triangles,viz.overestimation of the stiffness and sensitivity to distorted meshes.This strain smoothing approach has been widely employed to manifold engineering fields,for instance,nonlinear elasticity [20],coupled with XFEM [21] and hybrid/enrichment [20,22].A recent application to polytopes by Francis et al.[23,24] shows that strain-smoothed approach offers accurate and high performance in solution compared to both standard FEM and S-FEM.The salient features of S-FEM are: (a) does not require the gradient of the shape functions to compute the bilinear/linear form and (b) avoids Jacobian and iso-parametric mapping which lead to tolerate severe mesh distortion.Further,this scheme is relatively easy to implement to the existing FE codes with little modifications.In this study,the SIMP method is utilized within the framework of the cell-based linear smoothed finite element method for the sake of its computational reliability and robustness.The gradients are smoothed over the smoothing domains with a linear polynomial weight function,which is then employed to compute the stiffness.Materials are allocated on each smoothing domain and the weight function for the sensitivity filter is also defined to smoothing domain-based by the distance between each the smoothing domains.For this reason,the proposed strain smoothing topology optimization is fully smoothing domain dependent unlike the polyFEM.

This paper is organized as follows: Governing equations for linear solid elastic problems are briefly revisited in Section 2.Then,cell-based linear smoothed finite element approximation in the context of a polygonal finite element cell is introduced in Section 3.In the Section 4,the structural topology optimization within cell-based strain smoothing approach is explained.A few numerical examples are investigated in Section 5 to demonstrate the performance of the proposed scheme and followed by conclusion section.

2 Governing Equations

Consider an elastic body occupyingΩ∈Rdand it is bounded by (d－1)-dimensional surface.The boundary is assumed the following decomposition:Γ=Γu∪Γtand Ø=Γu∩Γt,where Dirichlet boundaryΓuand Neumann boundaryΓt.The boundary-value problem for linear elastostatic is: find u:Ω→Rd:

whereσis the Cauchy stress tensor,n is the unit outward normal f is the body force,and ^u and g are prescribed displacements and tractions on each Dirichlet boundaryΓuand Neumann boundaryΓt,respectively.

The following variational form is obtained from Eq.(1):

where the test function v.

The infinitesimal strainεis given as:

and its relation to the stress is given as follows using Hooke’s law:

where C is the elasticity tensor.

Substituting Eqs.(3) and (4),Eq.(2) can be rewritten as:

Then,the discrete form of Eq.(5) can be define by standard Bubnov-Galerkin procedure as follows:

where the stiffness K and the external force b are defined as:

where the shape functionsΨand the strain-displacement matrix from the derivatives of the shape functions B=∇sΨ,respectively.

3 Cell-Based Linear Smoothed Finite Element Approximation

In this section,the linear smoothing function in the framework of cell-based smoothed finite element method for polygonal finite element cells is briefly revisited.The infinitesimal strain tensor for smoothed strain approach on the smoothing domainΩkwith the smoothing functionf (x)is defined as:

In the standard strain smoothing approach,the smoothing functionfis a constantf (x)=1.In general,this scheme is used for linear triangular or bilinear quadrilateral finite element cells.However,the constant smoothing function is not suitable for polygonal meshes [23].Hence,an another treatment is required to use polygonal finite element in the framework of strainsmoothed scheme.Therefore,Francis et al.[23] proposed a linear smoothing function employing a polynomial basis function to the use of polygonal finite elements:

where the sampling point is given as x={x,y}andΨis the shape functions.Then,Eq.(10) can be rewritten as follows:

and its spatial derivative is obtained as:

where forΨI,1andΨI,2,respectively:

The aforementioned procedure leads to yield the basis function derivative level: therefore,the RHS of Eq.(8) can be rewritten as [23,24]:

whereΩkis the sub-domain,Γkis the boundary of sub-domain,n is the outward unit normal vector and the shape functionsΨfor arbitrary complex polytopes are the following Wachspress functions.

Wachspress basis function for polytopes.To generate shape functions of arbitrary polygons,Wachspress basis function [26] is used in this work.Wachspress introduced rational basis functions on arbitrary polygonal meshes that the algebraic equations of the edges are used to yield nodal interpolation and linearity on the boundaries.Later,Warren et al.[27] proposed the generalized Wachspress shape function for simplex convex polyhedral elements.For the interested readers,the detailed explanation of the use of Wachspress function for smoothed polygonal element cell can be found in [23,28].A simple form of a barycentric Wachspress basis function for irregular n-sided polygons is given by Meyer et al.[29]:

where

whereA(p,q,r)is the signed area of triangle(p,q,r),andγjandδjare adjacent angles as shown in Fig.1.Note that,the Wachspress basis function shows the lowest shape function while presenting the boundedness,linearity and linear consistency on convex polytopes.

Figure 1:Barycentric coordinates for a polygon: the circle is a field node and the triangle is a centroid

Modified basis functions in the framework of linear smoothed finite element approximation.The explicit form in two-dimension of Eq.(14) with the proposed the linear smoothing function can be explained as follows:

and

with the linear smoothing function derivativesf,1=[0,1,0] forΨa,1andf,2=[0,0,1] forΨa,2.

As defined in Eqs.(17) and (18),the proposed approach requires two different numerical integration schemes: one on the boundaries of smoothing domains (Γk) and another in smoothing domain (Ωk).For the integration on the boundaries,the aforementioned Wachspress shape functions are used,whilst the Lagrangian basis functions are employed for interior integration scheme.Two integration schemes are illustrated in Fig.2.As shown in Fig.2,the polygonal finite element cell is sub-divided into five triangular sub-cells with geometric centerO.As the same as the standard cell-based smoothed FEM,the target sub-cell is constructed as the smoothing domain in this approximation.Then,two Gaussian quadrature points are located on each boundaryΓkn,n={1,2,3} of the smoothing domainΩk.Three interior Gaussian quadrature points are in the smoothing domainΩkto compute RHS of Eqs.(17) and (18).The two-dimensional matrix form of the outward normal vector n is given as:

Now,Eqs.(17) and (18) can be rewritten as the two-dimensional system of equations:WGaussdj=fjwherej∈{1,2}.Gauss weights matrix WGaussand the vector fjare defined as:

and

where the coordinates ofm-th interior Gauss pointsmr=(mr1,mr2)and their weightsmw.The coordinates ofk-th boundary Gauss points and their weights are given asand,respectively.

Figure 2:Numerical integration scheme of the smoothing domain for the linear smoothing function for the polygon.A ‘filled’triangle is the geometric center of the polygon and ‘open’circles are field nodes.The interior Gauss points where the modified shape functions are obtained are depicted as red crosses and Gauss points on the boundary of the smoothing domain are given as‘filled’squares

Thej-th shape function derivatives vector dj,that is the modified shape functions and unknown vector in the present method,can be written as follows:

and djcan be obtained at three internal Gauss points as follows:

Then,the modified strain-displacement matrix can be evaluated as follows:

and the nodal strain-displacement matrix constructed atk-th internal Gauss point:

wherek∈{1,2,3}.

4 Solid Isotropic Materials with Penalization

In structural optimization,compliancecis minimized to attempt the optimal material distribution.To minimize the compliance,solid isotropic material with penalization (SIMP) method has gained massive attention among other topology optimization schemes,because of its use of fewer design variables and computation simplicity.

Therefore,the following power-law approach [30] is utilized in this paper to solve the problem in the context of the strain smoothing approximation:

where compliancecis twice the strain energy obtained by corresponding displacements,U and ueare the global and element displacement,respectively.NSDis the number of cell-based smoothing domain andis the local smoothed stiffness matrix built over smoothing domains.The penalization powerp=3.0 ～9.0 (in particular,the penalization powerp≥3.0 for Poisson’s ratioν=0.3),the volume fractionVfis the ratio of the material volumeV(x)and the design variable x is in this work is the density.Note that,the relative density xminis non-zero minimum value for the stability of simulation.The optimality criteria (OC) method [4] is employed to solve Eq.(22).The updated scheme using OC is given as follows [30,31]:

A tuning parameterηand a move limitζare used as control parameters to obtain stable convergence [31].,defined asBxeatk-th iteration,can be written as:

wherexerelies on the Lagrangian multiplierλdetermined by a bisection method.Eq.(22) can be written in terms of the design variables only:

wherec(xe)=fTu andThe compliance of Eq.(26) can be rewritten by adding the zero function:

In the current approach,the smoothed stiffness matrixis obtained over each smoothing domain.

Now,the sensitivity filter is discussed.The density-based topology optimization often faced from the notorious checkerboard patterns [32].Therefore,Bensøe et al.[31] introduced the following filter ensuring mesh-independency:

where the weight function/mesh-dependent convolution operatoris defined as:

Figure 3:The filtering scheme comparison between poly FEM and the proposed linear smoothed CS-FEM: (a) polyFEM and (b) polyCS-FEM.Closed triangles are the centroid of each polygonal cell and open circles are the centroid of cell-based smoothing domains

Figure 4:Topology optimization algorithm in the framework of the proposed linear smoothed cellbased strain smoothing scheme

5 Numerical Examples

In this section,a series of numerical tests are examined for the proposed topology optimization scheme.The domain is discretized with unstructured polygonal meshes using the open-source Matlab PolyMesher with builtin Matlab functions voronoin proposed by Talischi et al.[12].The current approach is compared to results of topology optimization with polygonal finite element cells introduced by Talischi and co-workers [10,13].A Matlab code for polygonal topology optimization written by Talischi et al.[13] is used as the reference solution to compare the performance of the current approach.The tuning parameterη=0.5 and move limitsζ=0.2 are often used for the rapid and stable convergence [31].The other parameters for the optimization in the following numerical examples are the same as those used in the reference [13].As the aforementioned,polygonal meshes are discretized to three-noded triangular sub-cells.Note that,the following examples are considered as dimensionless:

For discussing the results,the following convention is used:

· polyFEM: conventional finite element method with polygonal cells;

· polyCS-FEM: linear smoothing cell-based smoothed finite element method with polygonal cells;

· BCs: boundary conditions;

· DOFs: degrees-of-freedom.

5.1 Serpentine Beam

Firstly considering a Serpentine beam,for this test,the following control parameters are used:tuning parameterη=0.5,move limitζ=0.2,the volume fractionVf=0.5 and the penalization powerP=3.Fig.5 illustrates boundary conditions and a representative polygonal mesh generation which are used in this study.Young’s modulusE=105and Poisson’s ratioν=0.3 are used as material properties in this section.The external concentrated forceP=1000 is imposed as shown in Fig.5.Unstructured polygonal finite element are used with 10,136 DOFs.

Figure 5:Boundary conditions and unstructured polygonal cells for Serpentine beam: (a) red solid line is Dirichlet BCs (fully constraint) and red arrow is the external force,and (b) unstructured polygonal finite element cells

Obtained optimized topologies of the polyFEM and polyCS-FEM for the Serpentine beam are shown in Fig.6.The optimal topology of the proposed method is similar to its of Talischi et al.[13].The rate of strain energy reduction and the corresponding values are provided in Fig.7 and Table 1,respectively.For this test,the compliance reduce rate of the polyFEM is around 78%,while the polyCS-FEM is able to reduce only 77%,c=18.8572×103toc=419.84 with 86 iterations.

Figure 6:Optimized topologies of Serpentine beam: (a) polyFEM and (b) polyCS-FEM

Figure 7:The rate of compliance convergence for Serpentine beam

Table 1:Compliance (×103) at 1st and the last iteration for Serpentine beam

5.2 Suspension Triangle

A suspension triangle is studied in this example.The corresponding boundary conditions and mesh discretization (10,178 DOFs) are given in Fig.8.The same tuning parameter,move limit,volume fraction and penalization power as used in the previous example are utilized for this example.Isotropic material is assumed for Suspension triangle with Young’s modulusE=105and Poisson’s ratioν=0.3.The external forces are implemented at the left lower point as the concentrated forceP=1000 (see Fig.8).

Figure 8:Boundary conditions and unstructured polygonal cells for Suspension triangle: (a) constraint to x-direction at the left upper point,fully constraint at the left lower point and red arrows are the external forces,and (b) unstructured polygonal finite element cells

The optimized configurations for the suspension triangle derived from the polyFEM and polyCS-FEM are described in Fig.9.A notable performance of the proposed scheme for this test is obtained as shown in Fig.10 and Table 2 with the detailed values.After 165 iterations,value of 85.23% compliance reduction is observed by the current scheme.

Table 2:Compliance (×103) at 1st and the last iteration for Suspension beam

Figure 9:Optimized topologies of Suspension triangle: (a) polyFEM and (b) polyCS-FEM

Figure 10:The rate of compliance convergence for Suspension triangle

5.3 Multiple Load Problem

Next,multiple load problem is investigated in this section.Michell beam is considered in this section with two loadsP=1000,the corresponding boundary conditions and a schematic spatial discretization with polygonal meshes which is equivalent to 10,128 DOFs are shown in Fig.11.The following parameters are employed for the study: Tuning parameterη=0.4 and move limitζ=0.2.The volume fraction for this method is set toVf=0.5 and the penalization powerP=3.

Fig.12 represents the optimized topologies of Michell beam obtained by the polyFEM and the proposed cell-based linear S-FEM.As clearly shown in Fig.12,the proposed approach shows smoother optimal topologies while some areas of the polyFEM have rather rough surfaces.The compliance convergence curve for both methods is given in Fig.13 and the detailed compliance values are provided in Table 3.In the polyFEM,after 122 iterations,the strain energy is reduced fromc=3679.05 toc=855.45 which is 76.75% reduction.In case of the proposed method,it achieves 77.29% strain energy reduction with 201 iterations.At the same level of iteration (which is 122),the polyCS-FEM is able to decrease the strain energy toc=792.01.

Table 3:Compliance (×102) at 1st and the last iteration for Michell beam

Figure 11:Michell beam problem: (a) geometry and boundary conditions and (b) unstructured polygonal finite element meshes

Figure 12:Optimized topologies of Michell beam: (a) polyFEM and (b) polyCS-FEM

Figure 13:The rate of compliance convergence for Michell beam

5.4 Compliant Mechanism Synthesis

Lastly,a compliant mechanism synthesis problem is considered [33,34].The domain of force inverter problem is [0,3]×[0,1] and unstructured polygonal meshes are discretized with 10,120 DOFs as shown in Fig.14.Note that the given domain is only half of the whole domain since it is symmetrical.The following material properties are given as: Young’s modulusE=1.0 and Poisson’s ratioν=0.3.To minimize the negative horizontal output displacement,the input spring constantk1=0.1 and output spring constantk2=0.1 are utilized respectively.The volume is constrained to 30% with tuning parameterη=0.5 and move limitζ=0.2.The penalization factorp=3 is also used for this test.Fig.15 shows the optimal topologies of polyFEM and polyCSFEM.The convergence rate of the compliant mechanism design is given in Fig.16.As shown in Fig.16,the output displacements which are compliance in this example for each polyFEM and polyCS-FEM are obtained as: (a)c=0.0934～－0.9918 for polyFEM and (b)c=0.0927～－0.9187 for polyCS-FEM.

Figure 14:Optimized topologies of force inverter problem: (a) polyFEM and (b) polyCS-FEM

Figure 15:Optimized topologies of force inverter problem: (a) polyFEM and (b) polyCS-FEM

Figure 16:The rate of compliance convergence for force inverter problem

6 Conclusions

In this paper,the cell-based linear smoothed strain scheme for structural topology optimization and the method are presented.The proposed smoothed-strain framework has the following distinguishing features: (a) does not require the explicit form of shape functions and (b) two different boundary and interior integration schemes are employed.This leads to the standard CS-FEM can be utilized in polygonal finite element meshes.To use polygonal meshes in the context of smoothed finite element analysis,the polynomial basis function is employed as the smoothing function while it is a constant function in the standard S-FEM.The proposed polyCSFEM requires two Gauss points on the boundaries of each smoothing domain and three interior Gauss points.Wachspress basis functions are used the inside of the smoothing domain to evaluate the shape function of the present approach.The smoothing domain-based filter is successfully incorporated into an optimization algorithm.From a set of numerical tests,the following conclusions can be drawn:

· in general the polyCS-FEM needs more iterations to converge the compliance compared to the polyFEM because the bandwidth of the polyCS-FEM is wider than the polyFEM;

· the polyCS-FEM provides the relatively smoother optimal topologies compared to polyFEM;

· the polyCS-FEM is able to effectively yield compliance reduction rate even at the same level of iteration as the polyFEM.

To obtain the high resolution of optimal topologies,the employment of polytree meshes [35]will be considered as one of possible topics for the future communication.

Funding Statement: Changkye Lee,Seong-Hook Kee and Jurng-Jae Yee would like to thank the support by Basic Science Research Program through the National Research Foundation (NRF)funded by Korea Ministry of Education (No.2016R1A6A1A0312812).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.