Problem with solving the heat equation

1

Hello all,

This has been a problem that I have been trying to solve for a while but could not figure out what is the issue.

IBVP:
Cubic domain of 1mm x 1mm x 1mm
Initial Temperature: 300K
Dirichlet BCs: Temperature: 320K at boundary X=0, Y=0, Z =0
Neumann BCs: Heat flux of -1000 mW/mm^2 at boundary X=1, Y=1, Z=1
dt = 1e-5
NSteps = 1000

Using approximate material properties of steel.

I get the following output for temperature profile (Paraview)

HeatEquation
HeatEquation1131×781 132 KB

I also find that the temperature goes below the initial temperature, but this should not happen.

What could I be doing wrong when using CG1 elements? The following is the code:

#Using dolfin 2019.02.dev
from fenics import *
import time

t_tot           = 1.e-5
num_steps       = 1000
dt              = t_tot / num_steps
rho             =  8.e-6

Flux            = -1000.e0

Cv              =  0.5e6
T_init          = Constant(300.0)
Tref            = Constant(300.0)
K               = Constant(16.3e0)

X_Len           = 1.0
Y_Len           = 1.0
Z_Len           = 1.0

EX              = 8
EY              = 8
EZ              = 8

mesh            = UnitCubeMesh.create(EX, EY, EZ, CellType.Type.tetrahedron)

ScalarSpace     = FunctionSpace(mesh,'CG',1)

def X0(x, on_boundary):     return near(x[0],0.,DOLFIN_EPS) and on_boundary
def X1(x, on_boundary):     return near(x[0],X_Len,DOLFIN_EPS) and on_boundary

def Y0(x, on_boundary):     return near(x[1],0.,DOLFIN_EPS) and on_boundary
def Y1(x, on_boundary):     return near(x[1],Y_Len,DOLFIN_EPS) and on_boundary

def Z0(x, on_boundary):     return near(x[2],0.,DOLFIN_EPS) and on_boundary
def Z1(x, on_boundary):     return near(x[2],Z_Len,DOLFIN_EPS) and on_boundary

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)

AutoSubDomain(X0).mark(boundary_subdomains, 1)
AutoSubDomain(X1).mark(boundary_subdomains, 2)
AutoSubDomain(Y0).mark(boundary_subdomains, 3)
AutoSubDomain(Y1).mark(boundary_subdomains, 4)
AutoSubDomain(Z0).mark(boundary_subdomains, 5)
AutoSubDomain(Z1).mark(boundary_subdomains, 6)

dss     = ds(subdomain_data=boundary_subdomains)

bcTX0           = DirichletBC(ScalarSpace,      Constant(320.0),X0)
bcTY0           = DirichletBC(ScalarSpace,      Constant(320.0),Y0)
bcTZ0           = DirichletBC(ScalarSpace,      Constant(320.0),Z0)

bcT             = [bcTX0,bcTY0,bcTZ0]

def q(theta):
    return (-K*grad(theta))

T_n             = Function(ScalarSpace)
T_n             = project(T_init, ScalarSpace)

T               = TrialFunction(ScalarSpace)
vT              = TestFunction(ScalarSpace)

ThermalForm     = rho*Cv*(T-T_n)/dt*vT*dx  + Flux*vT*dss(2) + Flux*vT*dss(4) + Flux*vT*dss(6) - dot(q(T),grad(vT))*dx 
aT, LT          = lhs(ThermalForm), rhs(ThermalForm)

fileResults         = XDMFFile("Solution/HeatEquation.xdmf") # the output is written in a .xdmf file and could be open with Paraview
fileResults.parameters["flush_output"]          = True
fileResults.parameters["functions_share_mesh"]  = True

T = Function(ScalarSpace)
T = project(T_init,ScalarSpace)
T.rename("Temperature","label")

t = 0.

fileResults.write(T,t)

for n in range(num_steps):

    t += dt
    print('Step', n+1)
    solve(aT == LT, T, bcT)
    T.rename("Temperature","label")
    fileResults.write(T,t)
    T_n.vector()[:] = T.vector()


print("Solution Complete.")

Thank you!

2

Looks like your boundary conditions are not consistent with your initial conditions, so you’re effectively seeing Gibb’s phenomenon.

3

Hi Nate,

Could you tell me more about how I could correct this?
If boundary conditions, what should it satisfy?

Thank you!

4

see, e.g. here.