Topology Optimisation Heat Conduction Problem


71
views
0
19 days ago by
Otineb  

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

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.

Similar posts:
Search »