Hello everyone,

I’m having some trouble implementing the energy equation alongside with the NS, so I can create a simulation for **steady thermal flow** in some microchannels. Here’s my code:

So, this first part is just standard steady NS solving and everything works great.

from dolfin import *

import mshr

import matplotlib.pyplot as plt

from mshr import *mu = 1.05E-3

rho = 1000## Generate mesh with mshr

domain = (Rectangle(dolfin.Point(0.,0.), dolfin.Point(0.2E-2,1.2E-2)) +

Rectangle(dolfin.Point(0.2E-2,0.2E-2), dolfin.Point(1.8E-2,0.4E-2)) +

Rectangle(dolfin.Point(0.2E-2,0.6E-2), dolfin.Point(1.8E-2,0.8E-2)) +

Rectangle(dolfin.Point(0.2E-2,1.0E-2), dolfin.Point(1.8E-2,1.2E-2)) +

Rectangle(dolfin.Point(1.8E-2,0.2E-2), dolfin.Point(2.0E-2,1.4E-2)))

mesh = generate_mesh(domain, 100)

V = VectorElement(“P”, mesh.ufl_cell(), 2)

Q = FiniteElement(“P”, mesh.ufl_cell(), 1)

TH = MixedElement([V, Q])

W = FunctionSpace(mesh, TH)## Defining the boundary

inflow = ‘near(x[1], 0.0)’

outflow = ‘near(x[1], 1.4E-2)’

walls = ('near(x[0], 0.0E-2) || near(x[1], 0.2E-2) || near(x[0],2.0E-2) || ’ +

'near(x[1], 1.2E-2) || near(x[0], 0.2E-2) || near(x[0],1.8E-2) || ’ +

'near(x[1], 0.4E-2) || near(x[1], 0.6E-2) || near(x[1],0.8E-2) || ’ +

‘near(x[1], 1.0E-2)’)## Defining boundary conditions

bcu_inflow = DirichletBC(W.sub(0), Constant((0,0.1)), inflow)

bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow)

bcu_noslip = DirichletBC(W.sub(0), Constant((0,0)), walls)

bcs = [bcu_noslip, bcu_inflow, bcp_outflow]

#Define variational forms

v, q = TestFunctions(W)

w = Function(W)

u, p = split(w)## Equation

F1 = mu*inner(grad(u),grad(v))

dx + rhoinner(dot(grad(u),u),v)dxdiv(v)

- pdx - qdiv(u)*dx

J = derivative(F1, w)problem = NonlinearVariationalProblem(F1, w, bcs, J)

solver = NonlinearVariationalSolver(problem)

solver.parameters[“newton_solver”][“krylov_solver”][“maximum_iterations”] = 50

solver.solve()

So far so good. Next I tried to use the velocity *u* result from the previous simulation in a convective heat transfer case:

## Defined some constants

cp = Constant(4182) #J/kgK

k = Constant(0.6) #W/mK

s = Constant(1.0E5) #W/m3 source term## New Function space and T function

E = FunctionSpace(mesh, ‘CG’, 1)

T = TrialFunction(E)

r = TestFunction(E)

Here the boundary condition is only defined in the inflow as the other boundaries are the outflow and walls, which are dT/dx=0, a condition which I assumed to be naturally assumed with the variational formulation.

#boundary conditions

bct = DirichletBC(E, Constant(0), inflow)

Finally the equation:

F2 = rho

cpinner(dot(grad(T),u),r)dx - kinner(grad(T),grad®)*dx

L = inner(s,r)*dx

solve(F2 == L, T, bct)

The error message when I try this is:

**Expecting second argument to be a Function.**

I defined T as the Trial Function in this case, so I don’t understand the error. But I’m also second guessing my solution to the convective term: rho*cp*inner(dot(grad(T),u),r)*dx

I’ve been stuck in this for ages and have found no solution in the forum or the tutorial pdfs, as there is no specific case of steady thermal flow documented.

Thanks for reading!