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

aashish_kishore · June 20, 2024, 11:10pm

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!

hherlyng · June 21, 2024, 6:10am

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

aashish_kishore · June 21, 2024, 3:25pm

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

hherlyng · June 23, 2024, 9:28am

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

Topic		Replies	Views
Topology optimization for heat conduction using dolfin adjoint dolfin-adjoint	12	3503	November 11, 2023
Issue with a multidomain heat transfer variational formulation	9	323	March 25, 2024
Basic example does nothing variational formulation	1	441	January 26, 2021
Steady Thermal Flow with Navier Stokes, T is not a Function? variational formulation	7	683	January 29, 2021
Seperated Domains within problem General	3	514	August 1, 2023

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

Related topics