### Topology Optimisation Heat Conduction Problem

71

views

0

Hi,

I am trying to solve a topology optimisation problem for channelling heat. I have a 1x1 domain which I have equally divided into 3 sections on each side. From the bottom middle section (Dirichlet BC = 100deg C) I am trying to channel the heat to the top middle section (Dirichlet BC = 0deg C). The heat flux J is defined as dT/dn. However the code is not converging and material seems to coming from the outer sides of the middle section as shown below. The iteration was stopped at 100. Below is the objective function plot, solution and code. Does any one have any ideas why this is happening? Many Thanks.

```
from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *
import pyipopt
import numpy as np
tol = 1E-14
# turn off redundant output in parallel
parameters["std_out_all_processes"] = False
V = Constant(0.4) # volume bound on the control
p = Constant(5) # power used in the solid isotropic materil
eps = Constant(1.0e-3) # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8) # regularisation coefficient in functional
def k(a):
return eps + (1 - eps) * a**p # SIMP rule
n = 250
mesh = UnitSquareMesh(n, n)
A = FunctionSpace(mesh, "CG", 1) # function space for control
P = FunctionSpace(mesh, "CG", 1) # function space for solution
#Boundary Condition test 1
# ---------------------------------------------------------
class leftT(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],0,tol) and between(x[1],(0.6,1))
class leftB(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],0,tol) and between(x[1],(0,0.4))
class leftM(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],0,tol) and between(x[1],(0.4,0.6))
# ---------------------------------------------------------
class rightT(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1,tol) and between(x[1],(0.6,1))
class rightB(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1,tol) and between(x[1],(0,0.4))
class rightM(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1,tol) and between(x[1],(0.4,0.6))
# ---------------------------------------------------------
class topL(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],1,tol) and between(x[0],(0,0.4))
class topR(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],1,tol) and between(x[0],(0.6,1))
class topM(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],1,tol) and between(x[0],(0.4,0.6))
# ---------------------------------------------------------
class botL(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],0,tol) and between(x[0],(0,0.4))
class botR(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],0,tol) and between(x[0],(0.6,1))
class botM(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1],0,tol) and between(x[0],(0.4,0.6))
# ---------------------------------------------------------
boundary_markers = FacetFunction('size_t',mesh)
#leftM().mark(boundary_markers,1)
#rightM().mark(boundary_markers,1)
topM().mark(boundary_markers,1)
botM().mark(boundary_markers,2)
#bclb = DirichletBC(P, Constant(0), leftB())
#bclt = DirichletBC(P, Constant(0), leftT())
#bclm = DirichletBC(P, Constant(0), leftM())
#bcrb = DirichletBC(P, Constant(0), rightB())
#bcrt = DirichletBC(P, Constant(0), rightT())
#bcrm = DirichletBC(P, Constant(0), rightM())
bcbl = DirichletBC(P, Constant(0), botL())
bcbr = DirichletBC(P, Constant(0), botR())
bcbm = DirichletBC(P, Constant(100), botM())
#bctl = DirichletBC(P, Constant(0), topL())
#bctr = DirichletBC(P, Constant(0), topR())
bctm = DirichletBC(P, Constant(0), topM())
bc = [bcbl,bcbr,bcbm,bctm]
#bc = [bcbm,bcbl,bcbr,bcrm,bclm]
#bc = [bclb,bclt,bclm,bctl,bctr,bctm]
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
#dx = Measure('dx', domain=mesh, subdomain_data=domain_markers)
# ---------------------------------------------------------'''
#f = interpolate(Constant(1.0e-2), P) # source term
#g = Expression('-0.01', degree=1) # sample source term
#g = Expression('0.01', degree=1)
def forward(a):
T = Function(P, name="Temperature")
v = TestFunction(P)
F = inner(grad(v), k(a)*grad(T))*dx #- f*v*dx # with the addition of nueman boundary condition
solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
"maximum_iterations": 10}})
return T
if __name__ == "__main__":
a = interpolate(V, A) # name="Control") # initial guess.
T = forward(a) # solve the forward problem once.
controls = File("output/control_iterations.pvd")
a_viz = Function(A, name="ControlVisualisation")
J_cb = []
def eval_cb(j, a):
a_viz.assign(a)
controls << a_viz
J_cb.append(j)
np.savetxt('output_2d_l2/Objective_l2_v07.csv', J_cb, delimiter=',')
n = FacetNormal(mesh)
#J = Functional( f*T*dx + alpha * inner(grad(a), grad(a))*dx)
#J = assemble(dot(grad(T),n)*ds(1) - dot(grad(T),n)*ds(2) ) # objective function/ compliance
J = Functional( Constant(1000) * inner(grad(T),n)*ds(1) + inner(grad(T),n)*ds(2))
m = Control(a)
Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)
lb = 0.0 # lower bound
ub = 1.0 # upper bound
class VolumeConstraint(InequalityConstraint):
"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
def __init__(self, V):
self.V = float(V)
self.smass = assemble(TestFunction(A) * Constant(1) * dx)
self.tmpvec = Function(A)
def function(self, m):
reduced_functional_numpy.set_local(self.tmpvec, m)
# Compute the integral of the control over the domain
integral = self.smass.inner(self.tmpvec.vector()) # Compute the integral of the control over the domain
if MPI.rank(mpi_comm_world()) == 0:
print("Current control integral: ", integral)
return [self.V - integral]
def jacobian(self, m):
return [-self.smass]
def output_workspace(self):
return [0.0]
def length(self):
"""Return the number of components in the constraint vector (here, one)."""
return 1
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100}
solver = IPOPTSolver(problem, parameters=parameters)
a_opt = solver.solve()
xdmf_filename = XDMFFile(mpi_comm_world(), "output_2d_l2/material_v7.xdmf")
xdmf_filename.write(a_opt)
T_opt = forward(a_opt)
xdmf_filename2 = XDMFFile(mpi_comm_world(), "output_2d_l2/temp_v7.xdmf")
xdmf_filename2.write(T_opt)
```

Community: FEniCS Project

### 2 Answers

0

Hi,

May be stopped in a local minimum, try to use a continuation method in the SIMP relaxation start from 1 to 5 softly.

regards,

Francisco

May be stopped in a local minimum, try to use a continuation method in the SIMP relaxation start from 1 to 5 softly.

regards,

Francisco

0

Hi,

The problem of distributing discrete material values 0 1 in a fixed continuum domain to minimize the energy functional is an ill-posed problem, one way of solving this is by assuming a relaxed material model such as SIMP, as the penalizer (p) increases the behavior approaches of 0 1 design. I recommend you consult a more appropriate literature on this than this informal explanation.

I suggest earlier when possible, to know better its function in some aspects, analytical solution, if it has a global minimum, etc, from there to apply optimization techniques. In real problems it is difficult to obtain analytical solution beforehand.

You can consult references http://www.dolfin-adjoint.org/en/latest/documentation/poisson-topology/poisson-topology.html#bendsoe2003

http://orbit.dtu.dk/en/publications/numerical-instabilities-in-topology-optimization-a-survey-on-procedures-dealing-with-checkerboards-meshdependencies-and-local-minima(06b709ff-229c-49c6-89d3-05fad5fff501).html

regards,

Francisco

The problem of distributing discrete material values 0 1 in a fixed continuum domain to minimize the energy functional is an ill-posed problem, one way of solving this is by assuming a relaxed material model such as SIMP, as the penalizer (p) increases the behavior approaches of 0 1 design. I recommend you consult a more appropriate literature on this than this informal explanation.

I suggest earlier when possible, to know better its function in some aspects, analytical solution, if it has a global minimum, etc, from there to apply optimization techniques. In real problems it is difficult to obtain analytical solution beforehand.

You can consult references http://www.dolfin-adjoint.org/en/latest/documentation/poisson-topology/poisson-topology.html#bendsoe2003

http://orbit.dtu.dk/en/publications/numerical-instabilities-in-topology-optimization-a-survey-on-procedures-dealing-with-checkerboards-meshdependencies-and-local-minima(06b709ff-229c-49c6-89d3-05fad5fff501).html

regards,

Francisco

Please login to add an answer/comment or follow this question.

Thanks for the reply, do you mean change the parameter p? Could you explain?