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)
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!