Convection-diffusion equation: temperature plot has no signs of convection

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:
actual temperature plot

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!

Hi Aashish,

May it be that the scales of your square problem are such that convection does not play a dominant part in transport?

Please provide a MWE (Minimum Working Example), and I can try to help. The easiest would be to remove the Navier-Stokes part of the problem and just setting a velocity vector of e.g. (1, 0) to begin with, this would simplify the code a whole lot, whilst still letting us debug the code.

Cheers,
Halvor

Hello Halvor,
Thanks for your response!
I agree that scaling of the problem could be an issue, at the moment I have imposed a Reynolds number and then derived the velocity and pressure scales from Re.
Here is the MWE without the Navier-Stokes bit:

from dolfin import *
from ufl import tanh
tol = DOLFIN_EPS
mu = 1.004e-3
densityfluid = 998.0
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)
T_h = FiniteElement("CG", mesh.ufl_cell(), 2)
FlowSpace1 = FunctionSpace(mesh, U1_h)
TemperatureSpace = FunctionSpace(mesh, T_h)
TemperatureSpaceAdj = FunctionSpace(mesh, T_h)
u1 = Function(FlowSpace1)
t_fwd = Function(TemperatureSpace) 
t_adj = Function(TemperatureSpaceAdj)

# 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)
Pe_fluid = Constant(700)
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)


bct_walltop = DirichletBC(TemperatureSpace, Constant(350), Wall_top())
bct_wallbottom = DirichletBC(TemperatureSpace, Constant(300), Wall_bottom())
bc_ENE = [bct_walltop, bct_wallbottom]

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

ENE_fwd = derivative(ENE, t_adj, TestFunction(TemperatureSpace))
u1_out = File("u1_/plot_u1.pvd")
t_out = File("t_/plot_t.pvd")

assign(u1, interpolate(Constant((1.0, 0.0)), FlowSpace1))

# 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 << u1
t_out << t_fwd

Regards,
Aashish

Thanks Aashish! I think the reason you are seeing no convective behavior is because of the boundary conditions for the temperature. Dirichlet BCs are set on the top and bottom wall, and your variational form reflects zero Neumann BCs on the temperature on the inlet and outlet, so nothing else than the linear profile you’re observing is a possible solution. This is what I got by adding a non-zero Neumann BC on the outlet:

Here is your code with my edits:

from dolfin import *
from ufl import tanh
tol = DOLFIN_EPS
mu = 1.004e-3
densityfluid = 998.0
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)
T_h = FiniteElement("CG", mesh.ufl_cell(), 2)
FlowSpace1 = FunctionSpace(mesh, U1_h)
TemperatureSpace = FunctionSpace(mesh, T_h)
TemperatureSpaceAdj = FunctionSpace(mesh, T_h)
u1 = Function(FlowSpace1)
t_fwd = Function(TemperatureSpace) 
t_adj = Function(TemperatureSpaceAdj)

# 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)
Pe_fluid = Constant(700)
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)


bct_walltop = DirichletBC(TemperatureSpace, Constant(350), Wall_top())
bct_wallbottom = DirichletBC(TemperatureSpace, Constant(300), Wall_bottom())
bc_ENE = [bct_walltop, bct_wallbottom]

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0) # Create subdomains for the boundary integral measure
Inlet1().mark(subdomains, 1) # Mark the inlet by 1 (currently unused, does not appear in variational form -> flux is set to zero)
Outlet1().mark(subdomains, 2) # Mark the outlet by 2 
ds = Measure('ds', domain=mesh, subdomain_data=subdomains) # Boundary integral measure

ENE = inner(dot((u1/U), nabla_grad(((t_fwd-temp_2)/delta_temp))), t_adj) * dx \
    + (1/Pe_fluid) * inner(grad(((t_fwd-temp_2)/delta_temp)), grad(t_adj)) * dx \
        + inner(Constant(10.0), t_adj) * ds(2) # Impose a Neumann BC on the outlet

ENE_fwd = derivative(ENE, t_adj, TestFunction(TemperatureSpace))
u1_out = File("u1_/plot_u1.pvd")
t_out  = File("t_/plot_t.pvd")

assign(u1, interpolate(Constant((3.0, 0.0)), FlowSpace1))

# 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-6
solver_ENE_fwd.solve()

u1_out << u1
t_out  << t_fwd

With this change there is some convection of the temperature. Let me know if anything is unclear:)

Best,
Halvor

1 Like