Steady Thermal Flow with Navier Stokes, T is not a Function?

variational formulation
1

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)dx
- pdiv(v)dx - 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 = rhocpinner(dot(grad(T),u),r)dx - kinner(grad(T),grad(r))*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: rhocpinner(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!

2

The second input here should be a dolfin.Function, Which will contain the solution once the problem is solved.
I would suggest the following:

Th=Function(V)
solve(F2==L, Th, bct)
2 Likes
3

omg, you solved it! Thank you so much, you just saved me.

4

Apologies for butting in here, but this flow problem is similar to a dissolution problem I’m working on, so I was curious to get it running. The complete code (with correct markup and a couple of cosmetic mods) is this:

from dolfin import (Point, VectorElement, FiniteElement, MixedElement, Function, FunctionSpace,
                    DirichletBC, Constant, TestFunction, TestFunctions, TrialFunction,
                    split, inner, grad, dot, div, derivative, dx, solve,
                    NonlinearVariationalProblem, NonlinearVariationalSolver)
import matplotlib.pyplot as plt
from mshr import Rectangle, generate_mesh

mu = 1.05E-3
rho = 1000


# Generate mesh with mshr
domain = (Rectangle(Point(0.,0.), Point(0.2E-2,1.2E-2)) +
          Rectangle(Point(0.2E-2,0.2E-2), Point(1.8E-2,0.4E-2)) +
          Rectangle(Point(0.2E-2,0.6E-2), Point(1.8E-2,0.8E-2)) +
          Rectangle(Point(0.2E-2,1.0E-2), Point(1.8E-2,1.2E-2)) +
          Rectangle(Point(1.8E-2,0.2E-2), 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)


# Define 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+rho*inner(dot(grad(u),u),v)*dx-p*div(v)*dx-q*div(u)*dx
J = derivative(F1, w)


# Solve
problem = NonlinearVariationalProblem(F1, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['krylov_solver']['maximum_iterations'] = 50
solver.solve()


# Define 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 trial function
E = FunctionSpace(mesh, 'CG', 1)
T = TrialFunction(E)
r = TestFunction(E)


# Boundary conditions
bct = DirichletBC(E, Constant(0), inflow)


# Solve
F2 = rho *cp*inner(dot(grad(T),u),r)*dx-k*inner(grad(T),grad(r))*dx
L  = inner(s,r)*dx
Th = Function(V)
solve(F2==L, Th, bct)




#

But I get another error:

---> 79 Th   = Function(V)
    ...
    TypeError: Expected a FunctionSpace or a Function as argument 1

This is probably because of a naive mistake on my part, but I’d appreciate a tip to fix it.

5

As implied by the error message, V is not a FunctionSpace. In your code, you have defined V as a finite element.
If you would like to use a function space based on V,
I would make V_=W.sub(0).collapse()

1 Like
6

That was a fast reply! Thanks.

Well, I get further errors, so I’m beginning to think that when Pedro said the problem was solved, he meant in a limited sense.

7

Well to clarify, the solution proposed by dokken, while the comment itself was correct and helped me, it did not solve the error as the FunctionSpace variable was another. It should be:

Th = Function(E)

V is a vector space, and T is not a vector. So I advise that you use E.

8

Perfect! That worked. Thanks.