Hi! I’m using fenics 2017.2.0 (from conda) and I’m trying to match my results for a convective flow to the results from Christon et al. 2002 (https://www.researchgate.net/publication/247405116_Computational_Predictability_of_Time-dependent_Natural_Convection_Flows_in_Enclosures_Including_a_Benchmark_Solution). It’s a cavity problem where the left side is hot, the right side is cold, and the top and bottom are insulating. However, I’ve noticed that my equation for temperature is evolving in a purely diffusive manner, and I can’t figure out why i’m not seeing the dot(grad(T),u) term have any effect on my results. I believe this is what’s keeping me from matching the results from the paper, and I don’t really know how to deal with it. I’ve run this code out to 15000+ iterations, at which point eventually things blow up. I’ve implemented this a few different ways, but to make it easier to look through I figured I’d write everything in the form of the IPCS scheme used in the NS equation tutorial. Here’s my code:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
T =15 # final time
num_steps =15000 # number of time steps
dt = T / num_steps # time step size
# Create mesh and define function spaces
mesh = RectangleMesh(Point(0,0),Point(1,8), 40, 200)
V = VectorFunctionSpace(mesh, 'P', 2) #velocity mesh
Q = FunctionSpace(mesh, 'P', 1) # pressure mesh
Y = FunctionSpace(mesh,'P',1) # temperature mesh
# Define boundaries
u_0 = Constant((0.,0.))
u_n = project(u_0,V)
class side(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[0] < DOLFIN_EPS or x[0] >(1.-DOLFIN_EPS))
leftWall = 'on_boundary && near(x[0],0)'
rightWall = 'on_boundary && near(x[0],1)'
walls = 'on_boundary'
t = 0
# Define boundary conditions
bcu_noslip = DirichletBC(V, Constant((0., 0.)), walls)
bct_leftWall = DirichletBC(Y,Constant(0.5),leftWall)
bct_rightWall = DirichletBC(Y,Constant(-0.5),rightWall)
bcu = [bcu_noslip]
bct = [bct_leftWall,bct_rightWall]
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(3)
side = side()
side.mark(sub_domains,5)
ds = Measure('ds',domain=mesh,subdomain_data=sub_domains)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
theta = TrialFunction(Y)
tau = TestFunction(Y)
# Define functions for solutions at previous and current time steps
# u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
theta_ = Function(Y)
j = Constant((0,1))
Pr = 0.71
Ra =3.4E5
C1 = Constant(sqrt(Pr/Ra))
C2 = Constant(1./sqrt(Ra*Pr))
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
k = Constant(dt)
t_D = Constant(0.0)
theta_n = interpolate(t_D, Y)
# theta_n1 = interpolate(u_D, Y)
u_n = interpolate(Constant((0.0,0.0)),V)
#using IPCS method
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*C1*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = inner((u - u_n) / k, v)*dx \
+ inner(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ inner(p_n*n, v)*ds - inner(C1*nabla_grad(U)*n, v)*ds \
- inner(j*(theta_n), v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = inner(nabla_grad(p), nabla_grad(q))*dx
L2 = inner(nabla_grad(p_n), nabla_grad(q))*dx - (1./k)*div(u_)*q*dx
# Define variational problem for step 3
a3 = inner(u, v)*dx
L3 = inner(u_, v)*dx - k*inner(nabla_grad(p_ - p_n), v)*dx
#variational problem for temperature update (step 4)
F4 = inner((theta - theta_n)/k,tau)*dx + inner(dot(grad(theta),u_),tau)*dx + C2*inner(grad(theta),grad(tau))*dx - C2*inner(nabla_grad(theta),n)*tau*ds(5)
a4 = lhs(F4)
L4 = rhs(F4)
xdmffile_u = XDMFFile('velocity.xdmf')
xdmffile_p = XDMFFile('pressure.xdmf')
xdmffile_t = XDMFFile('temperature.xdmf')
vtkfileU = File('heat_u.pvd')
vtkfileP = File('heat_p.pvd')
vtkfileT = File('heat_t.pvd')
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A3) for bc in bcu]
[bc.apply(A4) for bc in bct]
n1=0
# Time-stepping
t = 0
import csv
hmax = mesh.hmax()
# with open('average_vel_file.csv', mode='w') as file1:
# writer = csv.writer(file1, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
while n1 < num_steps:
# Update current time
# try:
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1,'bicgstab', 'hypre_amg')
# plot(u_)
# plt.show()
# Step 2: Pressure correction step
b2 = assemble(L2)
solve(A2, p_.vector(), b2,'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
b3 = assemble(L3)
[bc.apply(b3) for bc in bcu]
solve(A3, u_.vector(), b3,'cg', 'sor')
# Step 4: Temperature Step
b4 = assemble(L4)
[bc.apply(b4) for bc in bct]
solve(A4, theta_.vector(), b4,'bicgstab','hypre_amg')
t += dt
if n1%10==0:
# plot(u_)
# plt.show()
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
xdmffile_t.write(theta_, t)
if n1%200==1:
vtkfileU << (u_, float(40))
vtkfileP << (p_, float(40))
vtkfileT << (theta_, float(40))
# z1 = grad(u)
# plot(z1)
# plt.show()
# Update previous solution
# u_n1.assign(u_n)
# theta_n1.assign(theta_n)
u_n.assign(u_)
p_n.assign(p_)
theta_n.assign(theta_)
t_hot = Expression('0.5', degree=1, t=t)
bct_leftWall = DirichletBC(Y,t_hot,leftWall)
bct = [bct_leftWall,bct_rightWall]
# dt = dt*1.01
# area = assemble(Constant(1.0)*dx(mesh))
average = dot(u_n,u_n)*dx(mesh)
average1 = assemble(average)
print('iteration: ',n1)
print('average vel: ',np.sqrt(average1)/4)
maxu = u_n.vector().norm('linf')
print('cfl approx: ',maxu*dt/hmax)
n1+=1
# Hold plot
# plt.colorbar()
plot(theta_)
plt.show()
plt.clf()
plot(u_)
plt.show()