On the ersatz material approximation in level-set methods

Abstract

Abstract. The level set method has become widely used in shape optimization where it allows a popular implementation of the steepest descent method. Once coupled with a ersatz material approximation [Allaire et al., J. Comput. Phys. 194 (2004) 363-393], a single mesh is only used leading to very efficient and cheap numerical schemes in optimization of structures. However, it has some limitations and cannot be applied in every situation. This work aims at exploring such a limitation. We estimate the systematic error committed by using the ersatz material approximation and, on a model case, explain that they amplifies instabilities by a second order analysis of the objective function.Mathematics Subject Classification. 49Q10, 34A55, 49Q12.

The topological ligament in shape optimization: a connection with thin tubular inhomogeneities

“…(1.8), as we assume thenceforth for simplicity (see Remark 2.3 below about this point). This is the well-known ersatz material method in shape and topology optimization: see for instance [2,11,30,50] about the consistency of this approach.…”

A level set-based structural optimization code using FEniCS

“…We use the framework of the ersatz material, which is common in level set-based topology optimization of structures; see for instance [4,48]. It is convenient as it allows to work on a fixed domain D instead of the variable domain Ω, but also can create instability issues as pointed out in [15]. The idea of the ersatz material method is that the fixed domain D is filled with two homogeneous materials with different Hooke elasticity tensors A 0 and A 1 defined by…”

“…----------------------------------------------------------------------12 def _main():13 Lag,Nx,Ny,lx,ly,Load,Name,ds,bcd,mesh,phi_mat,Vvec=init(sys.argv[1])14 eps_er, E, nu = [0.001, 1.0, 0.3] # Elasticity parameters15 mu,lmbda = Constant(E/(2 * (1 + nu))),Constant(E * nu/((1+nu) * (1-2 * nu))) Create folder for saving files 17 rd = os.path.join(os.path.dirname(__file__),\ 18 Name +'/LagVol=' +str(np.int_(Lag))+'_Nx='+str(Nx))19 if not os.path.isdir(rd): os.makedirs(rd) Line search parameters21 beta0_init,ls,ls_max,gamma,gamma2 = [0.5,0,3,0.8,0.8] Stopping criterion parameters 25 ItMax,It,stop = [int(1.5 * Nx), 0, False] Cost functional and function space 27 J = np.zeros( ItMax ) 28 V = FunctionSpace(mesh, 'CG', 1) 29 VolUnit = project(Expression('1.0',degree=2),V) # to compute volume 30 U = [0] * len(Load) # initialize U Get vertices coordinates 32 gdim = mesh.geometry().dim() 33 dofsV = V.tabulate_dof_coordinates().reshape((-1, gdim)) 34 dofsVvec = Vvec.tabulate_dof_coordinates().reshape((-1, gdim)) 35 px,py = [(dofsV[:,0]/lx) * 2 * Nx, (dofsV[:,1]/ly) * 2 * Ny] 36 pxvec,pyvec = [(dofsVvec[:,0]/lx) * 2 * Nx, (dofsVvec[:,1]/ly) * 2 * Ny] 37 dofsV_max, dofsVvec_max =((Nx+1) * (Ny+1) + Nx * Ny) * np.array([1,2]) (px,py,phi,phi_mat,dofsV_max) Define Omega = {phi<0} 42 class Omega(SubDomain): 43 def inside(self, x, on_boundary): 44 return .0 <= x[0] <= lx and .0 <= x[1] <= ly and phi(x) < 0 45 domains = CellFunction("size_t", mesh) 46 dX = Measure('dx') Define solver to compute descent direction th 49 theta,xi = [TrialFunction(Vvec), TestFunction( Vvec)] 50 av = assemble((inner(grad(theta),grad(xi)) +0.1 * inner(theta,xi)) * dX\ 0e4 * (inner(dot(theta,n),dot(xi,n)) * (ds(0)+ds(1)+ds(2))) ) 52 solverav = LUSolver(av) 53 solverav.parameters['reuse_factorization'] = True 54 #----------MAIN LOOP ----------------------------------------------55 while It < ItMax and stop == False: 56 Update and tag Omega = {phi<0}, then solve elasticity system. 57 omega = Omega() 58 domains.set_all(0) 59 omega.mark(domains, 1) 60 dx = Measure('dx')(subdomain_data = domains) 61 for k in range(0,len(Load)): 119 A = assemble( S1 * eps_er * dx(0) + S1 * dx(1) ) 120 b = assemble( inner(Expression(('0.0', '0.0'),degree=2) ,v) * ds(2)) 121 U = Function(V) 122 delta = PointSource(V.sub(1), Load, -1.0) 123 delta.apply(b) 124 for bc in bcd: bc.apply(A,b) 125 solver = LUSolver(A) 176 cx,cy = np.int_(np.rint([px[dof]/2,py[dof]/2])) 177 phi.vector()[dof] = phi_mat[cy,cx] 178 else: 179 cx,cy = np.int_(np.floor([px[dof]/2,py[dof]/2])) 180 phi.vector()[dof] = 0.25 * (phi_mat[cy,cx] + phi_mat[cy+1,cx]\ 181 + phi_mat[cy,cx+1] + phi_mat[cy+1,cx+1]) --------------------------------------------------------------------184 if __name__ == '__main__': 185 _main()…”

Topology Optimization Methods for 3D Structural Problems: A Comparative Study

“…The geometry in a fixed mesh is typically mapped to a mechanical model using either immersed boundary techniques (e.g., via X-FEM [111]) or density-based mapping (e.g., via relaxed Heaviside functions [2]) [112]. When using the second mapping, the stiffness coefficients are expressed in terms of the level-set function via approximated Heaviside functions, i.e., the ersatz material approximation is used [113,114], where the void material is replaced with a soft material. In addition, the LSF is commonly updated using a Hamilton-Jacobi (HJ) equation [2,3,47], 7 thus requiring the solution of a pseudo-time PDE equation.…”