### Topology Optimisation Heat Conduction Problem

71
views
0
19 days ago by

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 *
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)

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)
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

0
19 days ago by
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
Hi Francisco,
Thanks for the reply, do you mean change the parameter p? Could you explain?
written 19 days ago by Otineb
0
19 days ago by
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.