Hi everyone,
I am currently working on replicating a topology optimization problem using the method outlined in this paper I am employing the adjoint method and MMA within a 2D domain. However, in my resulting temperature plots, the contours seem to reflect only diffusive heat transfer, lacking the expected effects of convection and flow velocity.
Below, I have provided an MWE that solves the forward Navier-Stokes and energy equations from my optimization scheme for a single fluid. The boundary conditions are as follows: pressure inlet on the left, outlets on the right, no-slip on the walls, and Dirichlet temperature conditions on the top (hot) and bottom (cold) walls.
I anticipated a temperature plot similar to the one below:
However, my temperature plot from the MWE looks like this:
The temperature profile in the middle (along y axis) is linear, which should be parabolic. The velocity and pressure plots are non-zero and physical.
Here is the MWE:
from dolfin import *
from ufl import tanh
tol = DOLFIN_EPS
mu = 1.004e-3
densityfluid = 998.0
rhofluid = Constant(densityfluid)
lx = 0.0005
ly = lx
N = 50
mesh = Mesh(RectangleMesh(MPI.comm_world, Point(0.0, 0.0), Point(lx,ly), N, N, 'crossed'))
U1_h = VectorElement("CG", mesh.ufl_cell(), 2)
P1_h = FiniteElement("CG", mesh.ufl_cell(), 1)
T_h = FiniteElement("CG", mesh.ufl_cell(), 2)
A = FiniteElement("DG", mesh.ufl_cell(), 0)
FlowSpace1 = FunctionSpace(mesh, U1_h*P1_h)
FlowSpaceAdj1 = FunctionSpace(mesh, U1_h*P1_h)
TemperatureSpace = FunctionSpace(mesh, T_h)
TemperatureSpaceAdj = FunctionSpace(mesh, T_h)
DensitySpace = FunctionSpace(mesh, A)
w1_fwd = Function(FlowSpace1)
(u1, p1) = split(w1_fwd)
w1_adj = Function(FlowSpaceAdj1)
(v1, q1) = split(w1_adj)
t_fwd = Function(TemperatureSpace)
t_adj = Function(TemperatureSpaceAdj)
rho = Function(DensitySpace)
# Scaling parameters
L_scale = ly
Re_val = 100.0
Re = Constant(Re_val)
U_scale = (Re_val*mu)/(densityfluid*L_scale)
U = Constant(U_scale)
P = Constant(densityfluid*U_scale**2)
Pr_fluid = 7.0
Pr_solid = 3.5
Pe_fluid = Constant(Pr_fluid * Re_val)
temp_1 = Constant(350)
temp_2 = Constant(300)
delta_temp = 50
# Geometry
class Inlet1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0)
class Outlet1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], lx)
class Wall_top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], ly)
class Wall_bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
bcu1_wall_top = DirichletBC(FlowSpace1.sub(0), Constant((0.0, 0.0)), Wall_top())
bcu1_wall_bottom = DirichletBC(FlowSpace1.sub(0), Constant((0.0, 0.0)), Wall_bottom())
bcp1_outlet = DirichletBC(FlowSpace1.sub(1), Constant(0.0), Outlet1())
bcu1y_inlet = DirichletBC(FlowSpace1.sub(0).sub(1), Constant(0.0), Inlet1())
bcu1y_outlet = DirichletBC(FlowSpace1.sub(0).sub(1), Constant(0.0), Outlet1())
bcp1_inlet = DirichletBC(FlowSpace1.sub(1), Constant(U_scale**2*densityfluid), Inlet1())
bc_NS1 = [bcu1_wall_top,bcu1_wall_bottom, bcp1_outlet, bcp1_inlet, bcu1y_inlet, bcu1y_outlet]
# BCs declaration for Energy: Forward problem
bct_walltop = DirichletBC(TemperatureSpace, Constant(350), Wall_top())
bct_wallbottom = DirichletBC(TemperatureSpace, Constant(300), Wall_bottom())
bc_ENE = [bct_walltop, bct_wallbottom]
# Problem formulation
NSR1 = inner(dot((u1/U), nabla_grad(u1/U)), v1) * dx\
+ (1/Re) * inner(grad(u1/U), grad(v1)) * dx\
+ inner(grad(p1/P), v1) * dx + inner(div(u1/U), q1) * dx
ENE = inner(dot((u1/U), nabla_grad(((t_fwd-temp_2)/delta_temp))), t_adj) * dx\
+ (1/Constant(Pe_fluid)) * inner(grad(((t_fwd-temp_2)/delta_temp)), grad(t_adj)) * dx
UFL_State = NSR1 + ENE
NSR1_fwd = derivative(UFL_State, w1_adj, TestFunction(FlowSpace1))
ENE_fwd = derivative(UFL_State, t_adj, TestFunction(TemperatureSpace))
u1_out = File("u1_/plot_u1.pvd")
p1_out = File("p1_/plot_p1.pvd")
t_out = File("t_/plot_t.pvd")
assign(rho, interpolate(Constant(1.0), DensitySpace))
###### SOLVE NSR1_fwd
Jac_NSR1_fwd = derivative(NSR1_fwd, w1_fwd) # Jacobian
problem_NSR1_fwd = NonlinearVariationalProblem(NSR1_fwd, w1_fwd, bc_NS1, Jac_NSR1_fwd)
solver_NSR1_fwd = NonlinearVariationalSolver(problem_NSR1_fwd)
solver_NSR1_fwd.parameters['nonlinear_solver'] = 'snes'
solver_NSR1_fwd.parameters['snes_solver']['linear_solver'] = 'lu'
solver_NSR1_fwd.solve()
# Solve ENE
Jac_ENE_fwd = derivative(ENE_fwd, t_fwd) # Jacobian Energy equation
problem_ENE_fwd = NonlinearVariationalProblem(ENE_fwd, t_fwd, bc_ENE, Jac_ENE_fwd)
solver_ENE_fwd = NonlinearVariationalSolver(problem_ENE_fwd)
solver_ENE_fwd.parameters['nonlinear_solver'] = 'snes'
solver_ENE_fwd.parameters["snes_solver"]["report"] = True
solver_ENE_fwd.parameters['snes_solver']['linear_solver'] = 'lu'
solver_ENE_fwd.parameters['snes_solver']['relative_tolerance'] = 1.0e-4
solver_ENE_fwd.solve()
u1_out << w1_fwd.sub(0)
p1_out << w1_fwd.sub(1)
t_out << t_fwd
The docker image I am using is: Quay
I am unable to identify the error that is causing the diminishing of the convective effects. Any assistance would be greatly appreciated.
Thank you!
