Oscillations Heat Transfer Simulations

[Quentin]

Hi Everbody,

I am facing a small issue while simulating a simple heat transfer conduction problem of a bar with a Dirichlet condition at one end an d a Neumann at the other end .

As long I do not include the density and capacity close to the time derivative everything is Ok ( only with a term \frac{\partial T}{\partial t} instead of \rho c_{p} \frac{\partial T}{\partial t}), but when I introduce it I get

1. some oscillations close to the Dirichlet boundary condition ( even though this can be cured with quadratic elements )
2. a temperature field that does not evolve in time

see attached pictures

Here is the code I used

# density (kg.m^-3) 
rho = 8200.0
# capacity (J/Kg.K)
cp = 435.0
RhoCp = Constant(rho*cp)
# conductivity (W/m.K)
k = 11.5

# Boundary/Initial condtions -------------------------------------------------- 
# build plate temperature (K)
Tplate = 353.0
T_plate = Constant(Tplate)
# initial temperature (K)
T0 = Constant(293.0)
# imposed flux
Q = 100.0

# New version with better time stepping
T = 1.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size

l = 20.0
b = 4.0
h = 2.0



mesh = BoxMesh(Point(0, 0, 0), Point(l, b, h), 40, 8, 4)
tol=1E-8

class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < tol

class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return  x[0] > l - tol
    
def Bottom(x, on_boundary):
    return near(x[0], 0.)

def Top(x, on_boundary):
    return near(x[0], l)
    
# Define subdomains for boundary conditions
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryTop = BoundaryTop()
boundaryBottom.mark(sub_domains, 1)
boundaryTop.mark(sub_domains, 2)
    
facets = MeshFunction("size_t", mesh, 2)
AutoSubDomain(Bottom).mark(facets, 1)
AutoSubDomain(Top).mark(facets, 2)
ds = Measure("ds", subdomain_data=facets)
    
# Define variational problem
V = FunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

# define initial temperature
# 1) Constant one does the job
Tinit = interpolate(T0, V)
Tinit.rename("Temperature previous", "temperature profile at previous time step")

#Tplate = interpolate(T_plate, V)
BC = []
BC1 = DirichletBC(V, T_plate, boundaryBottom)
BC.append(BC1)


Neumann = Q*v*ds(2)

# Simplified one 
# 1) Without rho*cp
F = u*v*dx + dt*dot(k*grad(u), grad(v))*dx - (Tinit + dt*f)*v*dx -dt*Neumann
a, L = lhs(F), rhs(F)

# 2) With rho*cp
F = RhoCp*u*v*dx + dt*dot(k*grad(u), grad(v))*dx - (RhoCp*Tinit + dt*f)*v*dx -dt*Neumann
a, L = lhs(F), rhs(F)


vtkfile = File('VTK/Temperature.pvd')
u = Function(V, name='Temperature')
t = 0
for n in range(num_steps):
    print("current time is " + str(t) + " s")
    print("Completion is " + str(100*(t/T))+ " %")
    t += dt
    # Compute solution
    solve(a == L, u, BC)
    # Save to file and 
    vtkfile << (u, t)
    # Update previous solution
    Tinit.assign(u)
    # Update current time

I used the variational formulation described in

What is the variational form for a transient heat equation with Dirichlet, Neumann, and Robin boundary conditions?

Thank you Jorgen, I understand what you mean with the Neumann and Robin boundaries. I used a notation that was confusing, my apologies. I revisited the equations and corrected my mistakes - part of which was omitting density and specific heat capacity. I posted an update where the correct equation is provided. Thanks for the reply. Regrads, Willem

which seems correct to me .

I also got inspired from

Temperature jump near boundry condition

Dear colleagues, I am using the heat equation based upon the FENICS heat_eq tutorial as per following code: from __future__ import print_function from fenics import * import time T = 1.0 # final time num_steps = 1000 # number of time steps dt = T / num_steps # time step size # Import mesh and define function space mesh = UnitCubeMesh(12,12,12); V = FunctionSpace(mesh, 'P', 1) # Define boundary condition def left_boundary(x, on_boundary): if (x[1]<DOLFIN_EPS): retu…

but interestingly these constant are not accounted there

Any help/advice trick/would be very appreciated.

thanks,

Quentin

[Quentin]

I found out that I was just a unit problem. when all distance are in mm your density should be in kg/mm^-3 which introduces a 10^-9.

This topic could be closed I think